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On certain iterative methods for solving linear systems 
By 
A. S. HOUSEHOLDER* and F. L. BAUER** 


Among the iterative methods, or methods of successive approximation, for 
solving large linear algebraic systems’ of equations, those that have received the 
most attention recently have been the JAcoBI and the Gauss-SEIDEL, and their 
accelerated forms. A recent paper by GASTINEL, however, deals with a method 
of another class, which might be called a method of projection, and which in- 
cludes the method of steepest descent. In the present note an attempt will be 
made to describe such methods in general, and to give rates of convergence 
where possible. 


For solving the real or complex system 
(1) Ax=h, 


the matrix A being assumed nonsingular and of order n, let x, represent any 
iterate, and let 


(2) S,=% — %,, 1,=h—Ax,=As,, 


represent the residual and the remainder, respectively, where x itself is the 
true solution. By a method of projection will be meant one which requires at 
each step the resolution of s, into two components, one of which is required to 
lie in some subspace which might be selected only at the time, while the other 
component is s,,, and is required to be smaller than s, in some norm. Selection 
of a subspace is made by selecting a matrix Y,, whose columns are linearly 
independent and form a basis for the subspace. In practice Y, is generally of 
rank 1, hence a single vector y,, although methods of block relaxation can be 
interpreted as the use of a matrix of some higher rank g,<m. The methods of 
Gauss-SEIDEL and of relaxation can, in fact, be interpreted as particular methods 
of projection. 

Analytically, a method of projection is one for selecting a vector 4, of di- 
mension g,, or a scalar yw, if e,=1, such that 


(3) S,41 = 5, — Y, 4, 


is smaller than s, in norm. However, the simplest norms are ellipsoidal norms, 
and only those will be considered here. By an ellipsoidal norm of a vector s is 
meant the nonnegative square root of 


(4) I|s||*=s%Gs, 
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where G is any positive definite matrix. Since G is positive definite it is ex- 
pressible as a product 

(5) G= P" Pp, 

where P is some nonsingular matrix. Hence the ellipsoidal norm is simply the 
Euclidean length of Ps. If ||s||; is taken to represent the ordinary Euclidean, 
or spherical norm of s, then [3] 


(6) IslsSHPAllsiisi, — Isl] SUP lls Ilsils 
where ||P||; is the spherical or spectral norm of P, which is the positive square 
root of the largest proper value of G. 
The vector u, is therefore to be selected so that 
| s, — y u, || 
is minimized. But 
(s, — Y, u,)" G(s, — Y,u,) = s# Gs, — s#GY, (V#GY,)7¥# Gs, + wi ¥*GY,u,, 
where 
w, = 4, — (VEGY,)AY4Gs,. 


Hence the minimum is obtained when w,=0, hence when 


(7) u, = (Y,"GY,)7 YAGs,, 
and in this event 
(8) I s.||? — |] s,42||* = s# GY, u,. 


Evidently s, is not known, although 7, is. Hence to make the process feasible 
it is necessary that the positive definite matrix G, which will be fixed throughout, 
and the matrix Y,, be chosen so that Y,"G is the product of some matrix by A. 
If A itself is positive definite, G=A is an obvious choice. This will be considered 
later. Otherwise two choices suggest themselves. The first is to take 


(9) G= A" A; 
the second is 
(10) G=I, Y,=A*V,, 


which means that Y, is selected only indirectly through V,. 
The first choice, (9), gives 
u, = (YH AM AY,)7Y# Ay, 
11 
my Isll*—[lssaall? = 4 Yu, 
The method is essentially equivalent to the application of the method of pro- 


jection to the equivalent system 
A4Ax=A*"h, 


and needs no further discussion at this point. 

The second choice, (10) leads to 
(12) 4, = (VFA A*®y,)> ViK9,, 
(13) I] s,|]® — || Spal]? == 277 V, »,, 
where the norm is the ordinary Euclidean norm. 
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If V, is a single vector v,, the simplest choice is to make each v, equal to one 
of the unit vectors e;. These can be taken in rotation in accordance with the 
GauB-Seidel method. Otherwise, at any step the particular ¢; can be chosen 
that maximizes the right member of (12), or, more simply, the one that maxi- 
mizes |¢/7,|. This corresponds to the method of relaxation. 


In case the ¢; are taken in rotation, no lower bound for the right member 
of (13) can be given. But in case ¢; is chosen to maximize |e} 7,|, then 


Jer] = [Ill llr] & 2-9]]47 I] [N51 
Also 
ef A AM e;= ||A*%e,||*s || Al|?. 
Let 
(14) y =y(A) =||Al| ||4+|| 


be the spectral condition number of A [4]. Then 

IIs-11? — Ils,+a1]? = Ils,]l*/( 4), 
hence 
(15) ls-+all/Il 5.1]? S Bete (ny*)>. 


GASTINEL’s method [J] requires more arithmetic operations but no search. 
It is to take v, so that |v,|=e, and 


vr,= II*-lle- 


In the real case this means that each element of v, is +1 and agrees in sign 
with the corresponding element of 7,. But 


IIrolle = IIrol] = AUF Isl 
vy A AM, = |A*»,|/?s n || A||?, 


which leads to the same bound (45). 
The method of steepest descent takes 


whereas 


v,=7,. 
Then the right member of (13) can be written 
| || s.]]# (sf.A% A s,)?/[s9 s, sf (A* A)?s,], 
and if 
(16) y(A) =expe, 
then according to an inequality of KANTOROVICH [5], this cannot be less than 
|| s,||* sech* «. Hence 
(17) || S,42]|/||s,]| S tanh e. 
The norm is the ordinary Euclidean norm. As for the vectors 7,, by applying 
(2), (3), (10), and (42), it can be seen that 


Il742117/ I] 71]? ns lI 7.1 |? || 4 A” y,||? ||A%7,||-* me Dy 
5* 
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and an application of the Kantorovich inequality gives 


(18) I7,+|]/||,]| S sinh e. 
Hence the norms of the 7, need not decrease monotonically. On the other hand, 


from (17) it follows that 
IIs.|I/] sol] = tanh’, 


and since 

lI7-I/Ilroll S7Ils.I1/II soll 
therefore 
(19) l7-ll S7[|70l] tanh’ e. 

In case A is positive definite, take 
G=A. 

Then 
(20) u, = (Y"AY,)+ YF 1,, 
(24) IIs]? — [sya]? = 79" ¥, u,, 


where the norm is the ellipsoidal norm defined by 
(22) IIs]? =s"As. 


In the method of ‘‘block relaxation’’, the matrix Y, is made up of columns 
of the identity J, hence Y,¥ AY, is a submatrix of A, presumably with known 
inverse. If Y, is a single vector y,, the GauB-Seidel method takes these to be 
the ¢; in rotation. Consider, however, the method of relaxation, where each y, 
is the particular e¢; that maximizes |e/7,|. The argument is almost identically 
the same as before, both for this method and for GASTINEL’s, except that y* is 
replaced by y in (15). Hence for the method of relaxation and for GASTINEL’s 
method 
(23) II s-+all?/|]5-|]? S 4 — (my), 
with the further difference, however, that the norm is defined by (22). Also 
if the definition 
(24) y (A) = exp 2e 
of ¢ replaces (16), relation (17) appears for the method of steepest descent, again 
with the understanding that the norm is defined by (22). 

It cannot be guaranteed that in any individual step there will be a reduction 
of either the remainder 7, or the residual s, in the Euclidean, or spherical, norm, 
although by an argument analogous to the one that leads to (19) it can be shown 
that . 

(25) II s,[]s S|] Sol|sexp ¢ tanh’ e. 


Also the Kantorovich inequality can be applied to show that 


(26) l7-+2Ils/II7-lls S sinh e, 
with ¢ defined by (24). 
Simple geometric interpretations of the methods are readily available [/, 2]. 
Often the negative logarithm of the convergence ratio is used in the literature 
as a measure of the rate of convergence, since the number of steps required to 
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achieve a given degree of precision is roughly proportional to this. Applied to 
(15) and to (17), respectively, the results are, to the first order, (2m y*)“ and 2y~?, 
showing that the method of steepest descent converges at the worst approxi- 
mately 4m times as fast as the method of relaxation or that of GASTINEL at the 


worst. 
A generalization of the Kantorovich inequality, to be published soon in this 
journal, gives the following result: for y as in (16), if 


|v"7,| = [lel IIr|] sechn, 120 
then for the righthand member of (43) 
ler | 


vt AANy, = ||s||*sech*(e + m) 


and _ .erefore 


Wall Stanh(e+7). 


For the relaxation method and the method of GASTINEL, sechy=1 [yn or 
n= In(/n+ |/n—1). This gives, in place of (15), 


Mlsetall < y*(Ynt Yn—1)*—1 
IIsell 58 (n+ Yn—1)?+1 


and a corresponding result in place of (23). 
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Eine Fehlerabschatzung zum Ritzschen Verfahren 
fiir inhomogene Randwertaufgaben 


Von 


N. J. LEHMANN 


Die Aufgabenstellung entspricht einer ktirzlich in dieser Zeitschrift erschie- 
nenen Mitteilung von G. BERTRAM}. Hier sollen mit etwas moderneren Methoden 
fast schulmaBig bessere und vor allem verallgemeinerungsfahige Ergebnisse her- 
geleitet werden. 

Gleichzeitig wird die Gelegenheit benutzt und auf einige wichtige, aber kaum 
beachtete Zusammenhange zwischen Variationsproblem und Integral- bzw. 
Integrodifferentialgleichungen hingewiesen. 


1. Die Aufgabe und das Ergebnis 

Aus Griinden der Einfachheit wird wie bei G. BERTRAM das spezielle reell- 
wertige Problem 
(4) y’—g(*)y=r(x), y(0)=y(1)=0, g,rES, OSg(x%)SOQ 
vorgelegt. Dabei soll S alle auf <0, 1) stiickweise stetige Funktionen umfassen. 
Das Ritzsche Verfahren greift am zugehérigen Variationsproblem an: 
(2) J [y’*(s) +49(s) ¥°(s) +27(s) y(s)]ds>Min, yEZ. 
Integrationen sind immer iiber das Grundgebiet <0, 1) erstreckt; Z enthialt alle 
auf <0, 1> stetigen und etwa stiickweise stetig differenzierbaren Funktionen, die 
den Randbedingungen aus (1) geniigen. 


Man lést die Aufgabe (2) in einem Teilraum Z,€Z, dessen Elemente durch 
einen linearen Ansatz 


(3) Vn = 259. 9, €Z 


bestimmt sind. Als Koordinatenfunktionen werden hier wie in! die Eigenfunk- 
tionen der verwandten Differentialgleichung y’’+Ay=0, y€Z benutzt: 


(4) y, = \2sinvxx. 
Fiir den Fehler d, (x) der Ritzschen Naherungslésung des Variationsproblems 
(2) ergibt sich 
(5) |d,(x)| SE (x, n, Q) VJ redx, x€<0,i> 
mit 
E (x n Q) = 1 Q *(1—~*) ( ~ 2sin?»y 2x 
1 a, — aan a 


7m 1V3(n+1)? \4a, ™ IN heeyeore + *]-° (=e) 








1 BERTRAM, G.: Verscharfung einer Fehlerabschatzung zum Ritz-Galerkinschen 
Verfahren von Krytorr fiir Randwertaufgaben. Diese Z. 1, 135—141 (1959). 
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oder 
_ 1 [VQVx(1—-x) — 2sin*y 2x \t Q 
E, (x, n, Q) tk 7 (n+1)? + ae: yf ) ll V1+0Q/n2 +1]. 


Diese Schranken liefern an den Randern den Fehler Null und sind genauer und 
einfacher als die bisher bekannten (vgl. auch das Beispiel in Abschnitt 5). In 
E, (x, , Q) geht Q nur als Wurzel ein. Die beim Beweis benutzten iibersichtlichen 
Abschatzungsverfahren gestatten auch Differentiationseigenschaften der Koef- 
fizientenfunktionen in (1) voll auszunutzen. 


2. Integraldarstellungen fiir die Lésung der Variationsprobleme (2) in Z und Z,, 
und fiir die Fehlerfunktion 


Als Hilfsmittel wird die Greensche Funktion G(x, s) = 
Randwertproblem 


(6) M(y)=—y"=h(x), hES, yEeZz 
herangezogen. Sie lést (6) in der Form 
(7) y(x) =f G(x,s)h(s)ds 
und gestattet mit der zugehérigen, in Z positiv definiten Bilinearform 
(8) (uMv)=fu'v'dx, u,vEZ 
fiir jedes Element aus Z die Darstellung? 
(9) y(x) = (G(x,8)My), yez 
(der Pfeil weist auf die in der Form benutzte Variable). 

AuBerdem gehért G(x, s) fiir jede Variable zu Z. 

Setzt man nun einmal in der zum Extremalproblem (2) gehérigen Variations- 
gleichung 
(10) f[y’n’' + (gy +n)n]ds=(nMy)+JS(qyt+r)nds=0, y,n€Z 
n(s)=G(x,s)€Z(s) ein, so ergibt sich mit (9) die gesuchte Lésungsdarstellung 


x(1—s) xSs 
s(i—x)ssx 


(11) y (x) + f G(x, s) [9(s) y(s) +7(s)]¢s =08. 
Im Raum Z, leistet die Funktion 
(42) G,,(%, s) = y Hel) oe) €Z, fiir x unds 


die gleichen Dienste wie vorher G(x,s). Die g,(x) sind dabei nach (4) ortho- 
normierte Eigenfunktionen 

S PrGe4%*=5,0, SP PQ4x=(%,M y) =¥?774,,. 
Das und die Darstellung (3) der Funktionen y, liefert sogleich 


(43) n(X) =(Gu(*,8)M yn), Yn EZy- 


2? LEHMANN, N. J.: Eine Integraldarstellung fiir selbstadjungierte Randwertauf- 
gaben. Math. Nachr. 14, 129—156 (1955), insbes. Abschn. 4—8. Allerdings wird 
dort das allgemeine selbstadjungierte Problem M(y) =AN(y) behandelt. 

3 LEHMANN, N. J.: a. a. O., Abschn. 5. 
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Lést man jetzt (2) entsprechend dem Ritzschen Vorgehen in Z,, so bleibt 
(10) die Variationsgleichung — nur miissen die Lésung y= y,, und die Variation 
Elemente aus Z, sein. Damit ist 4(s)=G,(x,s) zulassig und ergibt iiber (13) 
analog zu (11) 


(14) n(x) + J G, (x, s) [9(s) ¥n(s) +17(s)]4s =0. 


Die Integraldarstellungen (11) und (14) fiir die exakte und die Naherungs- 
lésung des Problems gestatten sofort die Berechnung der ,,Fehler-‘‘Funktion 
d,=y—~y,- Man findet durch Subtraktion 


(15) d,,(x) + f G(x, s)q(s)d,(s)ds + f(x) =0 


mit f(x) =f [G(x,s) — G(x, s)] (¢(s) ¥n(s) +7(s)) ds. 


Diese Aufgabe ordnet sich vollstandig der bekannten Theorie polarer Integral- 
gleichungen unter‘. 


Danach gibt es eine Resolvente R(x, s; A), mit der (15) nach d, aufgelést wird 
(16) d,,(x) = f R(x, s; —1)9(s) f(s) ds — f(x). 


Zur Abschatzung dieser Funktion werden Zusammenhiange zwischen der 
Resolventen und der Greenschen Funktion herangezogen, die im nachsten Ab- 
schnitt bereitgestellt werden. 


3. Darstellung der Funktionen G(z, s) und R(x, 8; —1) 
mit Eigenfunktionen der Aufgabe? 


(17) y"+At(x)y=0 ye. 


Dabei geniigt es hier, ¢(x) € S, =O vorauszusetzen. Zu diesem Problem gehéren 
die Bilinearformen und Betrage 


(u,v),= ft(x)u(x)v(x)dx, |||, = )/(u, m),, u,vES 
(18) sowie 
(uMv) =f u'(x)v'(x)dx, INl||] = («aM u), u,veZ. 
Die erste Form ist dabei eventuell nur semidefinit, so daB es auch ,,Null- 
funktionen“‘ d€ S mit ||d||,—0, d@ ==0 geben kann. 


Trotzdem bleiben Dreiecks-, Schwarzsche und Besselsche Ungleichungen ver- 
fiigbar. Uberdies ist fiir jede solche Nullfunktion d und jedes y€ S 


(19) (4, y),=0. 


Die Eigenfunktionen y zu den Eigenwerten Ai, kénnen — soweit vorhanden — 
nach 


(20) (%, Vo)t = Oyo, (¥,M yt) _ 6,4, 





4 LEHMANN, N. J.: Bemerkungen zu einer Klasse polarer Integrodifferentialglei- 
chungen. Math. Nachr. 9, 45—50. — a. a. O. Abschn. 6 und 8. 
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orthonormiert verwendet werden. Alle J/ sind wegen 
i =(¥,M ¥5)/(¥%4, ¥5):>0 offensichtlich positiv. 


Damit besteht fiir G(x,s) und die zugehérige Resolvente R(x, s; A) die in 
beiden Variablen gleichzeitig absolut gleichmaBig konvergente Entwicklung 


(21) G(x,s)= > ae tei gaie + D,(x,s), R (4,8; 4) = Yad) + D,(%,s). 


D,(x, s) ist dabei ein "S6llthiedi* ohne Eigenwerte, der fiir jede Variable eine 
Nullfunktion ist. 

Diese Entwicklungen werden besonders fiir ¢(x)=q(x) und ¢(x)=1 bendtigt. 
Im letzten Fall ist y}(x)=9,(x), A}=»?2* und — wegen || y||?=f y*dx>0 bei 
yé€S, #0 — keine Nullfunktion vorhanden, so daB (21) in die bekannte Dar- 
stellung 


(22) G(x,8) = emt) Go(*) o(s) 


ary 
tibergeht. 
Diese Entwicklungen ermédglichen fiir die vorher benutzte Resolvente 
R(x, s; —1) sofort verschiedene einfache Abschatzungen. 


Da D,(x,s) eine Nullfunktion ist, folgt mit (19)|(20)|(21) und 44>0 


DAN} < yr LAlayl* 
(+1)? 4 (ag)? 


= f G*(x,s)q(s)\ds <Q f G*(x, s)ds = 2 + x*(4 — x)?, 








(23) ||R(x,s; —1)|P#=/S R2(x,s; — 1) q(s)ds= 2, 


d.h. 
(24a) R(x, — ts UP x( —%), *€<0,1). 


Insbesondere fiir groBe Q lassen sich leicht noch giinstigere Abschatzungen 
angeben. Hier soll dazu vereinfachend ¢(x)=0, *€<0, 1) nur fiir endlich viele 
Punkte zugelassen werden. Dann wird wegen (19) fiir Nullfunktionen d€ S 
fqddx=0, zufolge der Stetigkeitseigenschaften also d=0. In (21) verschwinden 
fiir dieses =g(x) auch die Zusatzterme D,(x, s), so daB 


viagra aph, a 
und mit (4f-+ 1)? Af >0 aus (23) 


(24b) |R(x,s; —1)||-s VG(«, x) = x(4 — x) 
erhalten wird. 








4. Abschatzungen fiir die Fehlerfunktion 
Aus der Integraldarstellung (16) erhalt man mit der Schwarzschen Ungleichung 
ZU (%,U)¢ 


(25) | 4,,(x)| S||R(x,s; — 1)[leli/ll, +17 (*)1- 
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Nach den Vorbereitungen im vorigen Abschnitt fehlen hierin nur noch Schranken 
fiir | f(x)| und ||/||,. Dabei ist nach (15) und (22) 


(26) f(x) =fG4(x,s)p(s)ds mit G4(x,s) =G(x,s) —G,(x,s) = > Sols} ols) 
und p=qy,+7rES. y 


Die Funktion y, ist die Lésung des Variationsproblems (2). 


Die Abschatzungsméglichkeiten fiir /(x) sind vielfaltig, hier wird lediglich 
ein bequemer Weg vorgefiihrt. 


Aus (26) folgt sogleich mit der Cauchyschen und der Besselschen Ungleichung 








_| J wls)Sopax — glx) tf <= , 
| #(x)| 4 JB }s[> er [ > Umeda 
(27a) <{ > BO} spas 


v=n+1 


und 


We=SofdxsOfPdx= State s) p(s) ds]*dx 


oo d 2 
=@ pa Usper s = ae Lf pdx)? s saci fPrdx 





vy=n+1 
d.h. 
(27b) Mlle aoe a Li Prax 
Dabei ist 
(28) [ferdx]*=[f (gy, +7)*dx]}tSVO[fqyidx}#+ [fredx}t 


grundsatzlich als bekannt anzusehen. Um eine lésungsunabhangige Fehlerab- 
schatzung zu gewinnen, muB noch {qy2dx ersetzt werden. Als Ausgangspunkt 
wird die zugehérige Variationsgleichung (10) mit y=y,, 7€Z, gewahlt; fiir 


7=Y%Y, War 
SyPdx+fqyndx=—Jry,dxs(frrdx)t(f yfdx)i. 


Mit der fiir y€Z bekannten Steklowschen Ungleichung f y’*dx>2°f y*d x folgt 
daraus 


(29) mf yndx+Sqyndxs[frrdx}*[f ypdx]t. 
Weglassen des zweiten Summanden links fiihrt auf 
(Syidx)is 4 (fr2dx)s 


und das in (29) mit fgy2dx<Qfy2dx zu 





fades 2 : fracas 2 fras. 
Q 
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In (28) eingesetzt ergibt sich damit schlieBlich 
0 2dxjis ——*¥__—. + 1] (f rd x) 
(30) (f p2dx) ¥ Vas + (f dx}. 
Aus den Gleichungen (24a), (24b), (27a), (27b) und (30) 14Bt sich sofort das 
SchluBresultat (5) ablesen. Man iiberlegt dazu, daB aus Stetigkeitsgriinden 


E,(x,, Q) auch ohne die in Abschnitt 3 gestellten Zusatzbedingungen fiir g(x) 
verwendet werden darf. 


5. SchluBbemerkungen und Beispiel 


Die vorgefiithrten Abschatzungen sind nur als beispielgebend zu werten und 
lassen sich an manchen Stellen verbessern. 


Zum Beispiel wurde bei der Einschrankung von |/(x)| und ||/||, auf S. 64 
besonders grob verfahren. Entscheidend ist dort der Summenwert: 


x y,pdx)? mit p=qy,+r. 


Aus der vethuniaiuanedialle ergibt sich 
(34) 2X (Sebdx) =f pdx— 2X (Sobax)? 


und hier sae der letzte Term sofort durch die Koeffizienten der Ritzschen 
Lésung y,= Lae, ausgedriickt werden. Gleichung (14) liefert mit den ent- 


sprechenden Bezeichmungen 


ya, 9,(*) =— 51 Pe (2) fepas dh. a,vat=— fg,pdx 
v=1 


| yp? m2 
und (31) damit sogleich 
4 4 
(32) [2 Ueda) = [Slay +1)8dx—at Datel 
Dieser Ausdruck tritt jetzt in Gleichung (5) an Stelle des Faktors 
1 Q 2d x\s 
——*—— + 1| (J dx 
Lae Vi+Q/a* - |u 
und. kann die Fehlerabschatzung wesentlich verbessern. 


Obwohl in (32) alle Terme prinzipiell bekannt sind, wird man fiir kleinere Q 
vereinfachend das Integral f p*d x nach (30) ersetzen: 


09 funteraen frees lie ahaa Hf [ee 


Weiter kénnen in (27) fiir p zusatzliche Differentiationseigenschaften voraus- 
gesetzt und dann die Konvergenz von >) »* f g,pdx ausgenutzt werden. Man 


gelangt so zu einer Schranke (5) mit E(x, n, 9)=0(-,). Andererseits kann 


q(x) auch als negativ zugelassen werden, bzw. es werden allgemeinere Probleme 
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M(y)= >)(—1)’(m,y)™=r(x) vorgelegt. Die entscheidenden Hilfsmittel der 
v=1 


Abschnitte 2—4 stehen in entsprechender Verallgemeinerung weiter zur Ver- 
fiigung (vgl. FuBnote 2). 
Zur Anwendung der Schranken nach (5) wird vor allem die Funktion 


S(x,n) =( > 2eintyn st 








4 
bendtigt : venta =” 
s S(s,1) S(s,2) S(x,3) | — S(x4) S(x,5) S(x6) |... | S(x,00) 
0 0 0 0 0 0 0 


0,1 0,268 0,170 0,113 0,0749 0,0491 0,0318 
0,2 0,374 0,165 0,0696 0,0463 0,0463 0,0401 
0,3 0,351 0,0992 | 0,0865 0,0691 0,0397 0,0323 
0,4 0,247 0,134 0,0976 0,0496 0,0496 0,0326 
0,5 0,171 0,171 0,0683 0,0683 0,0383 0,0383 


S(%, n) = S(1— %, n) 


























ooooco 


Oft geniigt aber schon eine Abschatzung fiir S(x, ), etwa 


S(xm)s(2 > 4) sate 3 mt4. 


vy=n+1 


Um einen Vergleich mit den von BERTRAM! gewonnenen Resultaten zu ermég- 
lichen, wird das gleiche Beispiel wie dort durchgerechnet 


y"’ —xy=e8—x*-—2, y(0)=—y(1)=0. 
Die Lésung ist y= x(1— x), die Ritzsche Naherung 
Yq = 0,258012sin x x — 0,000001 sin 2% x + 0,009 556sin 3 2 x — 0,000003 sin42 x. 
Mit E,(x, n,Q) =£E,(x, 4, 1) erhalt man als Fehlerschranken 


x,(1—x)=0 0,1 0,2 0,3 0,4 0,5 
ly as val SO 0,018 0,012 0,017 0,012 0,017. 


Nutzt man die bessere Abschatzung mit (32), (33): 


|4,(x)| <4 ee + S(x,n)| x 


V3(m +1)? 
x (8+ ppc Stee Le (f wena) 


fiir »=4 (man beachte y,= y2 sin yx), so folgt ohne wesentliche Mehrarbeit 
x,(1—x)=0 0,1 0,2 0,3 0,4 0,5 
|y—y| SO 0,0099 0,0065 0,0095 0,0071 0,0096. 
Der tatsachliche Fehler (bei x=0,1:0,0025) wird dabei nur noch um etwa den 
Faktor 4 iiberschatzt (bei BERTRAM noch um einen Faktor > 10). 


Dresden A 1 
Ernst-Thalmann-Str. 9 





i) 





(Eingegangen am 15. Oktober 1959) 
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Convergence of approximate eigenvectors 
in Jacobi methods 
By 
ROBERT L. CAUSEY* and PETER HENRICI** 


1. Introduction 


In this paper we are interested in studying the convergence of infinite products 
of unitary matrices connected with Jacobi methods [4] for computing eigenvalues 
of Hermitian matrices. We obtain bounds for the error in approximate eigen- 
vectors resulting from a finite number of steps of a Jacobi process. Also, we 
are able to answer certain questions concerning the representation of an arbitrary 
unitary matrix as an infinite product. Use is made of the classical results of 
Jacos! [5] and of the recent results of HENRICI [4] and Pore and Tompkins [8]. 
Basically, we adhere to the notation and terminology used in [4]. 

Consider the following computational algorithm. Let A=A,=(a,,) be a 
Hermitian matrix of order » with eigenvalues A, (r=1, 2,..., m). One calculates 
a sequence of matrices A,, Ay, ..., A4,=(a‘?),... which are unitarily similar to 
A by the recurrence relation 


(4.1) Ays, = UP A, U, (k =0,1, , af 


The U,=(ui*)) are special unitary matrices of order m. For every value of k 
there is specified a pair ,—=(t,,7,)=(t,7) of indices (we omit the subscript k 
in the sequel for notational simplicity) satisfying 1<1<jSm, such that the 2x2 
matrix 

(4 .2) V, = ( 


(k) (k) 
, 


k) k 
a 
which is a principal submatrix of U,, is unitary. All other elements of U, satisfy 


1, r=c 


(k) — oni 
(1.3) Uy 5, i rc, 


The matrices U, are completely determined by the pairs 2, and the 2 x2 unitary 
matrices V,. 





* This author’s work was begun at the Space Technology Laboratories, Los Angeles, 
under the sponsorship of the U.S. Air Force. It was concluded at the Lockheed 
Aircraft Corporation, Missiles and Space Division, Sunnyvale, California, under the 
Lockheed General Research Program. 

** This author’s work was done at the Space Technology Laboratories, Los Angeles, 
and was sponsored by the U. S. Air Force. 
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Any set of rules for choosing the U, is called a Jacobi method. A Jacobi method 
is said to be convergent if 


(1.4) A,>A (ko), 


for all A under the adopted set of rules, where A is a diagonal matrix whose 
diagonal elements are the A, in some order. Either of the quantities 


(1.5) $= late? 
T#C 

or 

(1.6) y= Max | al? 


may be used as a measure of the closeness of A, to A. For the Jacobi methods 
to be considered in this paper, a necessary and sufficient condition for convergence 
is that S,—>0 (k->0o) or, equivalently, 4,0 (k->00). See either [2] or [4] for 
the proof. If the infinite product of matrices 


(1.7) U= J] U,=U,U,U,... 
k=0 
converges, then 


(1.8) U*AU=A, 


and the columns of U are a complete set of normalized eigenvectors of A. We 
shall investigate the convergence of (1.7) for three convergent Jacobi methods 
defined below: (a) the classical Jacobi method, (b) the quasicyclic restricted 
Jacobi methods [4], and (c) the threshold cyclic Jacobi method [8]. Throughout 


the paper we shall restrict V, to be of the form 

cos #, — e'"sin d 

1. VY.=(_, ; rhe 
(1-9) , (masa cos 1 


By the classical Jacobi method we shall mean a Jacobi method with the 
following set of rules for determining the U, of (1.1). Choose 2, such that 


(1.10) | af? | =p, 
and choose the V, of (1.9) such that 
(1.11) alt+}) = 9, 
The equation (1.11) will always be satisfied if the relations 
(1.42) pa = arg alt, 
* 2| al * 


$s 1 


hold simultaneously. Since Jacobi dealt only with real symmetric matrices A 
and orthogonal matrices U,, the set of rules (1.10)—(1.13) constitute a generali- 
zation of his original method. 

A Jacobi method is called cyclic if in every segment of N=n(n —1)/2 con- 
secutive elements of the sequence {m,} every pair (p,g) (14 SP<qSm) occurs 
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exactly once. This is a special case of a guasicyclic Jacobi method with period K, 
in which the following condition is imposed. In every segment of K=>N conse- 
cutive elements of the sequence {m,} every pair (p,q) (1 p<qSn) occurs at 
least once. We next define a restricted Jacobi method. Let #f be an angle satis- 
fying (1.13) and belonging to the closed interval [— 2/4, 2/4], and let y>0 be 
a constant angle not depending on &. Then any method in which g, is chosen 
to satisfy (1.12) and #, is chosen such that 


(1.14) sign 0, = sign oF, | 9, | = Min (| dF |, y) 


is called a restricted Jacobi method with bound y. In [4] it was proved that any 
quasicyclic, restricted Jacobi method with suitable bound converges. Further- 
more it is shown in [4] that, if A has m distinct eigenvalues, any quasicyclic 
Jacobi method in which #,, p, satisfy (1.12), (4.13) and #f € [— 21/4, 2/4] converges 
quadratically provided the off-diagonal elements are already sufficiently small. 
More precisely, if 4nu,<d, then? 

(1.45) Prix S 2n(K —1) bsenxtormga yh, 


where | 
d= Min| A, — A,| 


The threshold cyclic Jacobi method is a modification of the cyclic Jacobi 
method due to Pope and TompKINs [8]. It depends on the choice of a sequence 
of positive threshold values 4, (v=0, 1, 2,...) such that ¢,,,<¢, holds for all ». 
A fixed cyclic ordering is adopted for the 2,; y, always satisfies (1.12), and for 
each y, one chooses 

* h) 
(1.16) =| Bf, for |apl=4 


0, for |a{P|<z,. 

When all | a‘)|<t#, (r+c), t, is replaced by ¢,,,. Pope and TomPKINs prove 
that the thresho!” rnethod converges if Jim t,=0. In the present paper we shall 
take the sequence {f,} to be a geometric progression; 7.e., we shall assume 
(4.17) t,=ag’ (0<q<1). 


2. Norms of Matrices 


In this section we summarize for convenience of reference certain known 
results concerning norms of matrices. Let H=(h,,) and G denote square matrices 
of order ». The norm (cf. [1]) of a square matrix H is defined to be a non-negative 
number %(H) satisfying the conditions 


(2.1) N(H)>0 for H+0, 
(2.2) N (cH) =|c|-N(H) for all scalars c, 
(2.3) N(H +G) <N(H) + NRG), 
(2.4) N(H G) <N(H) -N(G). 





1 Mr. E. R. Hansen has shown (oral communication) that the exponent in the 
exponential function may be replaced by /2nKd-yy,. 
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This definition is equivalent to that given by Ostrowski [7] for his ‘‘multi- 
plicative norm’’. By (2.4) 
N(H) = N(H 1) SR(H) - RUN), 
whence 
(2.5) N(I)S]1 


for any norm. There exist norms for which ¥(J)=1, an example being the 
familiar spectral norm 


(2.6) , (H) = Max (ee 


x* x 


Note that for any unitary matrix U, R,(U)=1. The Euclidean norm is defined 
by the relation 


= \3 
(2.7) M(H) =(D |e?) 
It is well known that, for all H +0, 


-4— &,(H) 
(2.8) wits oy S!- 


Inequalities such as (2.8) hold quite generally. In fact, it was shown in [7] that 
if N(H) and N’(H) are any two matrix norms, then 


(2.9) a(R, R’) < ae < b(M, X’) 


where a(X, MN’) and 5(R, N’) are positive constants which do not depend on H. 
It is easy to show that 
(2.10) Max |/,,| <2, (H), 


and by (2.9) there exists a positive constant b(9,, 3), independent of H, so that 
(2.11) | rye] SB (Ry, N) - RCH) 


holds for all 7 and c and for any norm. 


3. Preliminary Theorem on Infinite Products 


A sequence of xm matrices B,, B,,..., B,,,... with B,=(b") is said to 
be convergent if a limit exists for each element of B,,. That is, if lim bi") = 5, ., 
m-— > co 


then the matrix B=(b,,) is said to be the limit of the sequence {B,}. If the 
sequence {B,} is not convergent, it is said to be divergent. An infinite product 


(3.1) (I + £,) (J + Ey) (2 + Es)... 


2 In defining convergence of infinite products of matrices in general, it is necessary 
to exclude singular factors from the product and also to distinguish between singular 
and nonsingular limits. However, since we are dealing only with unitary matrices 
in the present paper, it is unnecessary to complicate the discussion with such dis- 
tinctions. 
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of Xn matrices is here* called convergent if the partial product 


(3.2) P, =I (I+E,) 


converges to a limit. We now prove the following 


Theorem 1. Let N(H) denote some norm of the matrix H. Let! the factors of 
the infinite product (3.1) be nxn unitary matrices, and let the series \N(E,) con- 
verge. Then (3.1) converges to a unitary matrix U. Furthermore, for all m, 


(3.3) N(U — P,) Ss 7, 2 (E,), 
where c(M)>0O depends only on N and n. 


Proof. We first establish that if the sequence {P,} has a limit U, then U 
must be a unitary matrix. Let P,—U+Q,,. We have that Q,,-0 (m->oo). 
Since each factor J+ E, is unitary, P,, is unitary for all m. Thus 


Px P,, = I = (U + Qu)* (U + Q,); 
or 
I — U*U=U*Q,,+ Q2U + OF On: 


The right hand side of the last expression tends to the null matrix as m-—>oo, 
and this shows that U is unitary. Let r>m and note that 


We shall prove that P, has a limit with the help of the following 

Lemma. // H,, Hz, ..., H, and G,, Gg, ..., G, are any two sets of n xn matrices, 
then the following identity holds 
(3.5) A, Hy... H, — G, Gq... Gy = (Hy — G,) Gy... Gp + 


+ H,(H, — G,)G,...G, +++» ++ Hy... Hy_s(Hy — G,). 


The truth of (3.5) is easily established by induction on p. The details are 
left to the reader. 


Letting H;=1+E,,,;, G;=I, and combining (3.4) and (3.5), we have 
(3.6) P— Py=PalEmsr+ (I+ Emit) Emte + 
+ (I+ Ents) (Z + Ente) Emits +++ + (LZ + Ems) --- (+ £,-1) E,]- 
Products of the matrices (J+ £,) are unitary and so, by (2.4) and (3.6), 
(3.7) Ny (PF, — Pu) SMa (Emer) +++ + (E,), 


since the spectral norm of a unitary matrix is 1. Now for an arbitrary norm %, 
there exist® a(M)>0 and b(R)>0 such that 


a(®) -R(A) SN, (A) S b(M) - R(A) 


3 The spectral norm &, is fixed; hence, the a(%,, R) and b(N,, R) of (2.9) do not 
depend on it. Therefore, for simplicity, we have here written a(®) for a(M,, R) and 
b(N) for b(R,, N). 

Numer. Math. Bd. 2 6 
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holds for all H. We have therefore that 


(3.8) a(R) - RF, — Py) SM (P, — Fn) 
and 
(3-9) X MW(Ey) SOM) DY R(E;). 
k=m+1 k=m+1 
Letting c(2)=5(N)/a(M) and combining (3.7)—(3.9), we find that 
(3.10) M(P— Py) <c(Ry Y R(Ey). 
k=m+1 


By hypothesis the series >) 9(£,) is convergent. Therefore, given any ¢>0, 
there exists an integer 4 >0 such that 


(3.14) 2 RE) < Fay cay ° 


p+l 
We see from (3.10) and (3.11) that 


(3.12) M (Pte — Fm) < Fey 
holds for all m= and for o=1, 2,.... But (3.12) is precisely the condition that 


(3.13) (Pate — Padre! <é 


holds for all y and c, by (2.11). Hence P,, tends to a limit U, as m-—>oo. Finally, 
letting 7->co in (3.10), we obtain (3.3) and this proves Theorem 1. 


4. Convergence of Eigenvectors and Error Estimates 


We first give a rather general result about convergence and estiraation of 
the error in all the eigenvectors. Following this we give estimates for the three 
Jacobi methods considered in this paper. 

Theorem 2. Let the nxn Hermitian matrix A have n distinct eigenvalues, 
and let the sequence of matrices (1.1) be generated by a Jacobi method such that 
the angles 3, satisfy (1.13), every ®,€ [— 2/4, 2/4], and such that the series 

co 

> St 

k=0 
converges, where S, is defined in (1.5). Then the infinite product (1.7) converges. 
Furthermore, for all m> mg, we have 


4 co 
(4.1) My (U — Pq) SF D St 
where 
(4.2) d=Min|4,—4,|, 


and where mz, is the smallest integer such that for all m>m, 


d 
(4.3) Sh< ‘ 
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Proof. We shall make use of the fact that, for some ordering of the eigen- 
values A, of A; 


(4.4) [4,—ap?| SS (r =1,2,...,m) 


hold for all k. The inequalities (4.4) may be derived from a theorem of LIDsKIi 
[6] or from an equivalent theorem of WIELANDT [9]. Alternatively, (4.4) follows 
directly from Theorem 1 of [4]. Note that, while (4.4) is valid for any value 
of k, the ordering of the eigenvalues A, is generally not the same for all k. How- 
ever, our use of (4.4) does not depend on the ordering of the A,. 


Let the U, of (1.1) be expressed in the form 
(4.5) U,=I+E,, 
and let E,=(el*)). By (4.2), (4.3) and (1.9) we see that the 2x2 principal sub- 


matrix of E, corresponding to V, is given by 


(4.6) eft) 2) _ (09s b,-—1 —e'sind, 
ef) elf) e~*%sin#, cos?,—1 }’ 
and that all other elements of E, are zero. From the definition (2.7) we get 
MN, (E,) Ss [2 (cos 0, — 1)? -+ 2sin? #,]!, 


and with the help of some standard trigonometric identities, we obtain 
| 6 





(4.7) MN, (E,) = 8! jsin = : 
By the definitions (1.5) and (4.6) 
1 
8 s — §. 
(4.8) Hp j2 } 
Hence, by (4.3), we have 
id d 
(4.9) Um < Va 4 <{ (m>m). 
Using (1.6), (4.13) and (4.2)—(4.4) it follows that 
2| af! 2Uk Hk 
(4.10) |tan 28,| = Ja — a] <3 =4 
rr 2 4 


holds for k>m,. Hence, by (4.9), |tan 20,|<1 and |20,|<2/4 (k>mp), since 
#,€[—2/4, 2/4]. For these small angles it is easily verified that 4 |tan(#/2)|< 
|tan 20|. Hence for k>m,, we have from (4.7), (4.8) and (4.10) 


i.) 2 gt gt 2 
(4.41) N, (E,) = 8! jsin 3 < = |tan 28,| << > aie d Sh. 


By hypothesis the series >) S} converges. Therefore ))2,(E,) converges and by 
Theorem 1 the infinite product (1.7) converges. Now applying (3.3) to the case 
where N=MN,, we see from (2.8) and (2.9) that 

(4.12) c(M,) = a(h,_®,) n), 

Combining (3.3), (4.11) and (4.12) we obtain (4.1) and this completes the proof 


of Theorem 2. 
6* 
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As before, we let N denote the number of elements on one side of the diagonal 
of an ” Xn matrix, viz. 


(4.13) N= jn(n —1), 
and for convenience we define 
(4.44) h=(1— NH, 


We let [x] denote the integral part of x, 7.e., the largest integer not exceeding x. 
Finally we let o=[N/q?] and we state our main result as follows. 

Theorem 3. Let the nxn Hermitian matrix A have n distinct eigenvalues, 
and let the sequence of matrices (1.1) be generated by either Method (a): the classical 
Jacobi method, Method (b): a quasicyclic restricted Jacobi method, or Method (c): 
a threshold cyclic Jacobt method in which the thresholds t, satisfy (1.17). Then 
the infinite product of matrices (1.7) converges. Furthermore, for Method (a) 
(4.15) Ny(U — Py) $2 (nS) e— — (m>m,), 
for Method (b) 

(4.16) R,(U — P,) < 4% 


where 0 is given by (4.23) below, and for Method (c) 
(4.17) Ry (U — Py) S ZF giro) (2mN yg) (m > m), 


' 
"Um (m>o), 





where in each case the integer my has the same meaning as it had in Theorem 2. 


Proof. It is easy to show that, for any step in a Jacobi method in which 
the parameters #,, p, satisfy (1.12) and (1.13), the following equation holds: 


(4.18) Shai = Sp— 2 af? |*. 
In the classical Jacobi method we have 
(4.19) |af/P2EN7S, 


always satisfied. Hence 
(4.20) SriiS(1i-N4)S,=HS, (k=0,1,2,...) 
holds for Method (a). It follows therefore that 





co 
4 Tae a on a 
(4.24) DISH+A+H+--) = Shy, 
whence by Theorem 2 the product (1.7) converges. By (4.20) 
a 4, am 
2, ‘ts % 1—A 


holds for all m. For m>m, we may combine the last expression with (4.1) and 


obtain the inequality (4.15). 
Now consider Method (b). We note that u,—>0 as k->oo, since a quasicyclic 
restricted Jacobi method is convergent. Since A has m distinct eigenvalues, 


there exists a smallest integer #, such that 
Min (|Of|,y)=|OF] (RZ A), 
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and by (1.14) the restricted method behaves like an ordinary quasicyclic method 
from the step p, on. There also exists a smallest integer ,=>, such that 4nu,<d 
(k=), and hence (1.15) holds for k=>,. Since the sequence {S}} is monotone 
decreasing and ju, —>0, there is a smallest integer #, such that (4.3) is satisfied, viz. 


St< 4 (kz), 


and therefore (4.11) holds fork=>,. Finally, a smallest integer p, exists such that 


(4.22) Mh<2e* (k2%), 
where g is the coefficient of u? in (1.15). Now let 
(4.23) @ = Max (hq, ps, 4). 


For k=o0 we have that (1. 7 and (4.11) hold, whence 


pi MN, (E,) << — >) My 


Hs: 


ae fd as 


i=0 j7=0 
Ss 5 tg + gud + etuh + ht 
+ (Mes + 8 Mos + Bo Moga +0) te 
+ (Mot K—-1 + 8 Mor K—1 + 8 Me+K-1 + ***)} 





4a. 2 & by 8 Mo+1 8 Mot K-1 
< “ eee ox 
al dg 1—g 4, . 1-8 Mo+1 ¥. * 1—§ 6+ K-1 : 
and using (4.22) we get 

2 
(4.24) x, (Ey < 4 2B (Mg + Meta + +++ + Me+K-1)- 


k=@ 
Since S, is monotone decreasing, 


Meir S75 Ste S75 SN eS 5 ff, (v=0,14,2,...), 


and we combine this with (4.24) to obtain 


(4.25) 3° (Ey) < = te: 

k=e 
Now, combining (3.3), (4.12) and (4.25), we get (4.16) and, since ,,—>0 (m-—>oo), 
this proves the convergence of (1.7) for Method (b). 

Turning our attention to Method (c), we see by (1.16) that U,=J, E,=0 
whenever |a‘)|<t,, so that a step in the process (1.1) counts only when a non- 
trivial transformation is made, i.e., whenever |a{})|=t,. We let the scalar « in 
(1.17) be equal to uy, so that 


14,-| SM@o=Mo P= (r+). 


There may be no non-trivial transformation made at any given threshold level 
but, if some are made, we can easily bound their number. Assume that after 





) 
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the r-th transformation u,<t,_,. Then 


(4.26) S,<2N#, 
and, by (4.18), 
(4.27) Sri1 SS, — 28 


holds for each transformation made at the »-th threshold level, since | a{|2>?,. 
Combining (4.26) and (4.27) we see that 


S,4pS2NG_,— 208 


holds after ~ transformations at the v-th level. Letting [x] denote the integral 
part of x, we have that there are at most 





ot 


steps at the »-th level, for 


S,+0, < 2N #-.— y a 





|é 





<2N@_,—2(“8=- —1)¢ 


<2#, 


so that u,,,,<#,; t.¢., the next threshold level has been reached. From (1.17) ’ 


«=| =o (v=1,2,...). 


We have proved that 
S,.22N@ (vy=1,2,...), 
so that, using (1.17), 


(4.28) Sh, S (2N py)”. 

Now, since s[t/s]S# for all positive integers s and ¢, we may write 
° oo co 

(4.29) x St s2 Shiney) 


k=m 
oo 


S XL Shue: 
k=a([m/o) 


Let 4 be one of the integers 0, 1, ...,0—1, and let 


k=o(™|+i+jo. 


l-b+lel+ 4 


We have 


and, since i<o, we have 


[=| =1+[4 (i=0,1,...,0—1). 
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Thus we may write (4.29) as follows: 


o- 


co 1 © 
2 St = a Shi tion)» 
k=m t=0 j7=0 

co 


=? 2 Shi tmio})> 


and using (4.28) we have 
¥ S}<o(2N pop)t S git ov 
k=m 7=0 
1 
1—q- 
Finally, combining (4.1) with (4.30), we get immediately the expression (4.17). 
Obviously, 


=o gi™'"\(2NV y) 





g@"1+0 (moo), 


and thus the product (1.7) converges also for Method (c). This completes the 
proof of Theorem 3. 


5. Remarks 


If we refer to unitary matrices of the general type used in the process (1.1) 
as elementary unitary matrices, then it is a rather easy matter to show that an 
arbitrary unitary matrix can be decomposed into a product of a finite number 
of elementary unitary matrices. In [3] (footnote 5, p. 8) GIVENS gave a proof 
that every orthogonal matrix with positive determinant can be expressed as a 
product of a finite number of elementary orthogonal matrices. We shall not 
give the details here, but a proof of the corresponding result on the decomposition 
of an arbitrary unitary matrix can be constructed by a method very much the 
same as that used in [3]. 

In spite of the result just mentioned, it is perhaps not without some interest 
that we can use Theorem 3 to prove that an arbitrary unitary matrix can be 
expressed as an infinite product of elementary unitary matrices. The decompo- 
sition is, of course, not unique. We sketch a proof as follows. 

Suppose W is any Xm unitary matrix. Let D be a fixed xn diagonal 
matrix with m distinct and real diagonal elements. Let the Hermitian matrix 


(5.4) A=WDW* 
_ be diagonalized by (say) Method (a) of Theorem 3, so that 
U*AU=D", 


where U is the limit of the infinite product (1.7). The matrix D’ has the same 
diagonal elements as D but possibly with a different arrangement. By a finite 
number of elementary unitary transformations with V, of the form 


0-1 
(to): 
we can transform D’ into D, so that 
(5.2) V*AV=D, 
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where V is an infinite product of elementary unitary matrices. From (5.1) and 
(5.2) we see that 


(5.3) VDV*=WDW*. 


It is easy to show from (5.3) that W can be expressed as V Y where Y is a diagonal 
unitary matrix, 7.e. 


(5.4) Y,,=e* . (y= 14,2,...,%). 


Since any diagonal matrix Y satisfying (5.4) can be decomposed into a product 
of at most m —1 elementary unitary matrices with V, of the form 


eF oO 
Oo. é8)’ 


Acknowledgement. The authors wish to express their thanks to Professor GEORGE 
E. ForsyTHeE for calling their attention to some of the problems discussed in this paper. 


the desired result is established. 
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Bemerkungen zur Fehlerabschatzung 
beim Ritz-Galerkin-Verfahren nach Krylow 


Von 


W. BORSCH-SUPAN 


G. BERTRAM [1] hat in dieser Zeitschrift kiirzlich eine Verscharfung der Fehler- 
abschatzung von KryYLOow [2, 3] zum Ritz-Galerkin-Verfahren bei gewéhnlichen 
Randwertaufgaben zweiter Ordnung gegeben. Auch die Bertramsche Abschatzung 
kann noch erheblich verbessert werden, insbesondere wenn man, was oft der 
Fall ist, tiber Koeffizienten und rechte Seite Differenzierbarkeitsvoraussetzungen 
machen kann. Eine solche Verbesserung gibt auch N. J. LEHMANN [4], die im 
folgenden beschriebene Vorgehensweise weicht jedoch davon ab?. 


Die Grundidee der Krylowschen Abschatzung, der auch BERTRAM folgt, 
besteht darin, zundchst einmal den Fehler /,(x) einer n-gliedrigen Ritzschen 
Naherung zu vergleichen mit dem Abbruchfehler F,(x) einer nach n Gliedern 
abgebrochenen Reihenentwicklung nach dem System der Ansatzfunktionen?. So 
erhalt man im Falle der Randwertaufgabe 


y’—q(x)y=r(x), y(0)=y(1)=0, 
wo 
O0Sq(x) SQ 


ist, mit einem Fourier-Sinus-Ansatz die Abschatzung 
|fn(*)| S| Fi, (*)| + 4,() - [Fall 


welche die ersten drei und den Anfang des vierten Abschatzungsschrittes bei 
BERTRAM zusammenfaBt. Zur Abkiirzung ist 





Hale) = (*0=" + 06(8)) oes 





(1) 04(*) = 2 S Sees. 


=n+1 





1 Der Verfasser erhielt von dieser Arbeit erst bei der Niederschrift seines Manu- 
skriptes Kenntnis. LEHMANNSs Ergebnisse sind daher hier nur vergleichend im Zahlen- 
beispiel am SchluB der Arbeit verwertet. LEHMANN hat die Differenzierbarkeits- 
eigenschaften dabei aber absichtlich noch nicht beriicksichtigt. 


2 Die Bezeichnungen der vorliegenden Arbeit sind dieselben wie in [7] und werden 
der Kiirze halber hier nicht mehr genauer erlautert. 








80 W. Borscu-Supan: 


1 

gesetzt. Dabei steht ||F,|| fir ( S\E, (x)? x), und » muB der Ungleichung 
(n+ 1)?2?> Q geniigen. ° 

Scharfere Abschatzungen dieses Types, bei denen H,,(x) mit wachsendem n 
sogar gegen Null strebt, wurden von KryLow und BocorjuBow [2, S. 21] ent- 
wickelt. Unter der Voraussetzung zweimaliger Differenzierbarkeit von g(x) erhalt 
man an Stelle von H,,(x) den folgenden gegeniiber [2] noch leicht abgewandelten 
Ausdruck: 


(2) H3 (x) =Q(1+ ) (on(x) + 278 —2)_) 4 








1—Q 2-*(n+1)-? V3 22 (m+)? 
Zila’ ll+4ll¢"ll Q *(1—+) 
+ Moe (tt Sy) 


der fiir n—>oo gegen Null geht. Auch sind entsprechende fiir x=0 und x=1 
verschwindende Ausdriicke leicht zu erhalten. 

Im zweiten Teil des Abschatzungsprozesses gilt es, Schranken fiir den Rest F, 
der Fourier-Reihe der Lésung y zu finden. Hierbei wird sowohl in [17] wie auch 
in [2, 3] Genauigkeit verschenkt, weil die in vielen praktisch wichtigen Fallen 
vorhandene Glatte der Funktionen g und 7 nicht ausgeniitzt wird. Diese Glatte 
auBert sich in einem mit wachsendem m rascheren Abklingen der Reste der 
Fourier-Reihe von 7 und damit von y”’ und y. 

Wir fiihren den besonders iibersichtlichen Fall genauer aus, in dem sowohl 
q als auch 7 zweimal differenzierbar sind und quadratisch integrierbare zweite 
Ableitungen besitzen. An geeigneten Stellen verwenden wir bereits in [1] formu- 
lierte Abschatzungen wie 


1 1 ’ 1 
(3) IvlSse lil. to's Wl 


los (1 +S) Ill 
und die Ungleichungen 


(4) I <e>all 5S ere aye ll <e" Dall S apepaye Hel 
und 
(5) |<8>n| Son (2) |I<e" Dall. 


die giiltig sind fiir Funktionen g(x) mit quadratisch integrabler zweiter Ab- 
leitung und g(0)=g(1)=0. Die Klammern < 5, bezeichnen wie in [1] den Ab- 
bruchfehler der n-ten Fourier-Partialsumme. 


Aus der Differentialgleichung gewinnen wir 


I<9">oll = Fell S1<99>all + 1<>ol- 


In beiden Termen der rechten Seite soll die Glatte der vorkommenden Funk- 
tionen ausgenutzt werden. Dazu verwenden wir (4). Im zweiten Term ist 7(x) 
zuvor aufzuspalten in eine am Rande verschwindende Funktion 7*(x) mit r*”’=7" 
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und eine die Randwerte enthaltende lineare Funktion 
r**(x) =7(0) + (7(1) —7(0)) x, 
deren Abbruchfehlernorm direkt berechnet wird: 


y=n+1 v=n+1 
» gerade y ungerade 


(6) Keel Fe (ro —ray S L4fro+eayye > *). 
Damit entsteht 
WES ag cape QD)" I+ ata I+ eal 
Hier differenzieren wir im ersten Term das Produkt aus, ersetzen g und 


seine Ableitungen durch die Maxima ihrer Betrage und schatzen die Normen 
von y und seinen Ableitungen gem48 (3) ab. Dann finden wir 


W'S Mce** Dal + array [le + (@ (8+ Se) + 2 + Se) el 


wobei zur Abkiirzung 
(7) O= pax fal=)}, C= max lel), "= gman |e"(a)) 
gesetzt ist. 


Zu Abschatzungen fiir |F,| und ||F,|| fihren uns jetzt die Ungleichungen 
(5) und (4). 

Fiigen wir die beiden Teile des Abschatzungsprozesses zusammen, so gelangen 
wir zur Endformel 


Hf (*) 1 
ee [all + saeETE 


x {Ill +(@ (4+ $)+22+4) lel]. 


|f.(*)| S Jo, (x) + 





in der man sich (1), (2), (6) und (7) eingesetzt zu denken hat. 


Zur Vereinfachung kann man o,(x) und ||<7**>|| durch ein wenig gréBere 
obere Schranken von durchsichtigerer analytischer Form ersetzen. Diese erhalt 
man beispielsweise, wenn man den Wert 1 an Stelle von sin*»zx einfiihrt und 
die verbleibenden Summen mit der Eulerschen Summenformel abschatzt. Auf 
diese Weise kommt man zu 


1 2 3 : 
(8) ou(2) Sse Vacca (+ amery + eee) 


(9) I|<r**>, || S + Vas (1 + = + ea) V[r(0)]? + [r(4)]?. 




















Das asymptotische Verhalten von /,(x) ist nun deutlich zu erkennen. In 
den beiden Faktoren der rechten Seite iiberwiegen die an erster Stelle stehenden 
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Summanden, so daB 


W. Borscu-SuPAan: 


[fn (*)| =O(n-*) 


ist. Eine starkere Abnahme von |/,| mit wachsendem m kommt nur bei ver- 
schwindenden Randwerten von r vor. 

AbschlieBend vergleichen wir unsere Formel mit den von BERTRAM [1] und 
LEHMANN [4] erzielten Resultaten an dem dort gegebenen Zahlenbeispiel. Wir 
stellen die Ergebnisse in Tabellenform zusammen. 
































Tabelle 1 
copie Abschatzungen fiir den Fehler nach 
, ‘a icher 
— Fehler LenMann [4]* 
fa (x) BERTRAM [1] Borscu-SuPaNn 
(5) (32) (33) 
1) 0 0,0521 0 0 0,0000 
0,05 0,0028 0,0524 0,0159 0,0089 0,0045 
0,075 0,0030 0,0526 0,0184 0,0102 0,0052 
0,1 0,0025 0,0527 0,0178 0,0099 0,0050 
0,2 —0,0007 0,0530 0,0116 0,0064 0,0031 
0,3 —0,0012 0,0533 0,0171 0,0095 0,0046 
0,4 0,0902 0,0535 0,0128 0,0071 0,0033 
0,5 0,0015 0,0535 0,0172 0,0096 0,0046 


Verwendet man die Ungleichungen (8), (9) und x(1— x) <4} 


die Abschatzung 


| f4(%)| S0,0057. 


Auf eine Tabellierung der Funktion o,(x) wird verzichtet, da sie bis auf 
einen Faktor mit der von LEHMANN [4] tabellierten Funktion S(x, n)=2°o,, (x) 
iibereinstimmt. Jedoch erfaBt diese Tabelle fiir gréBere m nicht das Maximum 
der Funktion. Wir erganzen daher die Tabelle in [4] nachstehend noch um einige 
Werte des Arguments x aus der Gegend des absoluten Maximums. 


, so erhalt man 


























Tabelle 2 
erm S(x,1) S(x,2) S(x,3) S(x,4) S(x,5) S(x,6) 
0,02 0,068 0,051 0,042 0,036 0,031 0,028 
0,04 0,128 0,093 0,073 0,060 0,050 0,042 
0,06 0,182 0,127 0,095 0,073 0,057 0,045 
0,08 0,228 0,152 0,108 0,078 0,056 0,040 
0,10 0,268 0,170 0,113 0,075 0,049 0,032 
0,12 0,302 0,180 0,111 0,067 0,040 0,026 
0,14 0,329 0,184 0,104 0,057 0,033 0,027 
Maximum 0,379 0,184 0,113 0,078 0,058 0,045 
Abszisse des 

Maximums 0,227 0,143 0,104 0,081 0,067 0,057 











3 Die Werte dieser Spalte kann man nur mit Hilfe der Ritzschen Naherung, also 
erst nach der Lésung der Aufgabe berechnen, wahrend alle anderen Abschatzungen 
a priori gemacht werden kénnen. 
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On the efficiency of certain quasi-random sequences 
of points in evaluating multi-dimensional integrals 
By 
J. H. HALTON 


1. Introduction 


In a recent paper on Monte Carlo Methods [J], J. M. HAMMERSLEY considers 
the estimation of the multi-dimensional integral 


1 1 
Jar. [deat ) (1) 


by formulae of the form 2; w;/(x,,,..., *;,). He points out that the efficiency 
of such an integration formula may be gauged by considering how it fares when 
f(x, ..-, %,) is the indicator-function of the hyperbrick defined by an arbitrary 
point A=(A,,..., A,) in the unit hypercube U*. The value of the integral (1) 
would then be the volume V=A,A,...A, of the corresponding hyper-brick. 


He uses the efficiency criterion 
1 1 
J = faA,...f44,|S(4) — NVI, (2) 


where N is the number of points used in the estimating formula and S(A) is 
the number of these points falling within the hyper-brick defined by A, and 
discusses the relative efficiencies of different distributions of points. 


K. F. Rotu [2] shows that a positive constant c, exists, such that 


J>c,(InN)* (3) 


for all distributions of the N points in U. 


J. C. VAN DER CorPuT [3] has shown that a sequence of N points can be 
constructed, for which 


where C is a constant and the supremum is taken over all points (A,, A.) in the 
unit square. 

HAMMERSLEY has proposed a generalisation of VAN DER CORPUT’S sequence 
to k dimensions, and asks whether sup | S(A) — N V| can be shown not to exceed 
some multiple of a power of In N. 





* The unit hypercube is defined by 0S 4%; 51 (i=1, 2,...,%). The hyper-brick 
is defined by 03%;S4; (t= Ride vase k). 
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In the present paper, it is shown that, for HAMMERSLEY’s sequence, both 


sup|S(A) — NV|<C,(In N)*-} (5) 
and 
J < B, (in N)?*~*, (6) 


where B, and C, are known positive constants. The sequence may be made 
rather more convenient for continuous computation at the expense of replacing 
k—1 by & in (5) and (6). 


2. The quasi-random sequences in k dimensions 


Let R be any integer; then any other integer m can be written in radix-R 
notation as 
N = Ny Ny_}.-. Ne My Mo = Ny +n, R+ ny R2+.---+ ny R™, (7) 
where 
M = [logr] = [Inn/In R}, (8) 


square brackets denoting the integral part, as is usual. By reversing the order 
of the digits in m, we can uniquely construct a fraction lying between 0 and 1, 


P = Ye(t) = 0° gM, Ng... My = My R14 n, R24 +--+ + ny, Ro! (9) 


VAN DER CorPUT’S sequence of points in the unit square is (m/N, p,(m)) for 
n=1,2,...,.N. For this, he shows that the result (4) holds. HAMMERSLEY’s 
suggestion is to use the k-dimensional sequence, 


(n/N, pr, (), Pr,(), +++» Pr, (")) for n=1,2,...,N, (10) 


in which he takes R,, R,,..., R,_, to be the first k—1 primes. Although this 
is clearly the most useful suggestion, we shall only require in what follows that 
these integers be prime to each other. 

In using a sequence such as (10) for the computation of integrals of the 
form (1), it is necessary to compute the function /(x,,..., x,) at each point 
of the sequence. Unfortunately, the first co-ordinate in (10) depends on N, so 


- that we should have to decide in advance on the value of N in order to perform 


the calculation. It is therefore suggested that a sequence of the form 
(Pr, ("), Pr,(™),---»Pr,()) for n=1,2,...,N, (11) 


_ might be more convenient in practice, at least in establishing the magnitude 


required for N in using (10). 


3. Some preliminary results 


Choose an arbitrary fraction A =0-@)4,4,... ay... lying between 0 and 1. 
By (9), if A>, one of the following conditions must be satisfied: 


Ag >My; Ag=M, AM; A=—M, A=—M, A> Me; ae (42) 
ag =, a, = 1,.---,4y=Ny; 


where we assume that A is given in non-terminating form. 
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Returning to the original integer n, we see that (12) is equivalent to requiring 
that, for some m< M+}, 


n=a,+a,R+---+a,,.R"-*+n,,_,R”* (mod R"), 
where (13) 
Any > My, and ny,,=0. 


It is clear that only one of the conditions (13) can hold for any given m and A. 
Further, for a given value of m, (13) means that » must be congruent to one of 
a,,—, distinct numbers (only one in the sole case of m=M-+ 2). Let p represent 
any one of these numbers. Then the number of integers in the sequence 1, 2, ..., N 
which are congruent to , modulo R™, is [N/R™]+-h, where h is a variable whose 
only possible values are 0 and 1. ° 

We may now extend these considerations to several variables. R is replaced 
by a set of mutually-prime integers R;. Each integer m gives rise to a set of 
fractions y;=pz,(m), and we wish to count the number of such sets satisfying 
the simultaneous conditions A;>q,;, where A is an arbitrary point in the cor- 
responding unit hypercube. We are thus concerned with the problem of satisfying 
sets of simultaneous congruences, 


n = p, (mod R¥*). (14) 


By the Chinese Remainder Theorem [4], since the moduli R7* are prime to each 
other, the solutions of (14) form a complete congruence-class, modulo /7,R%™. 
Thus the numbers of integers in the sequence 1, 2,..., N which satisfy (14), 
and so which generate points ~=(g,) lying in the hyper-brick defined by A 
(the meaning of all A;>g;), will be [N/I7;R%™]+4h. To get all such points, we 
must sum over the possible choices of the m; and the #;. 


4. The efficiency of the sequences 


We may now proceed to calculate J for the sequences. The case of (11) is 
the more straightforward. Remembering that there are a,,_, integers p for each 
dimension (except when m= M +2, where now M=[log,N]) we see that 


M+1 k 
S(4) = E (11 Pima) (NIT RE] + 0), (15) 


where O[ 0S 1* and }; ,,=4; ,, except for 5; 4,41=1, 4; being the radix-R; digits 
of A;. Now, noting that, for each 1, 


A; = 2 4; m—1 RY. (16) 
m=1 
we observe that 
co k 
NV = S(T 4im-1) (NUTR. (17) 
Each m=1 =1 


Thus the error is 


M;+2 k co k 
|S()— NV =| X (1D Pims}0— (IT ame evi, (18) 








* There will be a distinct value of @ for each set m,, mg, ..., My. 
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where curly brackets denote the fractional part, and since [N/R@+!L]=0 for 
any integer L. Since each of the sums above is positive, we increase the expression 
on the right of (18) by increasing the greater sum or decreasing the smaller. Thus * 








M; k 
|S(A)—NV|smax(| 3 (LT bim-1) (0 — (HTRTY’). 
Mab (19) 


po TT ¢1.m-1) (8 — {N/T R)) 


Each mg=1 \b=1 
where the # and {N/IT; R7*}’ vanish when any m,;=M,+-2 and c, ,,=4; ,, except 
for ¢; m,41=4i,m,4i+1. Therefore, since | — {N/II; R7"}|<1 and the same result 
holds for the primed quantities, and since ¢; ,,=4; », 20, 


Mi+2 k k .Mi+2 
[s(4)—NV]<_ > {ET ems) = IT ( ms): (20) 
whence 
k 
sup |S(A) — NV| < J] [M,(R;—1) + (2R;—1)] < 
0 Pee (24) 
< (n.wy* TT (Fe) if N>every R;. 
Again, by (2) and (20) **, 
k 1 M,+2 2 k M,+2 2 
J<II faa,( a Cy) = IT 2 a | >» csm-1) . (22) 
t=10 m=1 t=1 Die m=1 


The squared quantity is one plus the sum of M;+ 2 integers, each taking one of 
the values 0, 1, 2,...,(R;—1). Any particular sum gq will occur 7, times, where 
T, is the coefficient of 2 in (4 f2+22+.---4+2%i-1)Mit+2, We thus see that 
k te 
— -—e e d 
J<[[R;™ * 2,7 (q +1)? = [] Ri - la (z as 2,%24)| 


t=1 


. I] Rr “2 (ee edtete takyry)| 


k 
= [1 {t+ 5 G+ 2)(R—1) +4 +2) (M+) (R—9+ F (23) 


++ (M, + 2)(R,— 1) (R,— 2}, 


Final ‘ k . 
J<4*T] (Mi +3) (Ri —1)9< 4-4 TT SER (InN +3 in Ry 


t=1 








* If the first sum is the greater, | S(4) — N V| will be less than the result of taking 
the second sum in the expression on the right of (18) up toM;+ 1 only. If the second 
sum is the greater, we correspondingly take the first sum up to M;+ 1 only, and 
put c for a in the second sum and terminate it at M;+ 2 (This is simply the process 
of “rounding-up” the fractions A,). It follows that |S(A)—NV| will in either case 
be less than the greater of the two resulting expressions. The inequality (19) will 
follow if we make all sums up to M;+ 2 by using the primed quantities as shown. 


Ry-1 
** Here > = a and c;,, and a;,, are related as previously stated. 
Possible Each 4; ,_4=0 { ; 


Digits 
Numer. Math. Bd. 2 7 
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Thus, when N > max (R?), for instance, 


I<(2 “TT Se) any (24) 


We may now turn to the case of the Hammersley sequence (10). The process 
described for k dimensions above now holds for (k — 1) dimensions only. Of the 
N fractions n/N (n=1,2,..., .N), only the first NA fall below A, and so only 
the first [NA] integers can be considered as candidates for a (k — 1)-dimensional 
sequence of type (11). Noting that [[NA]/L]=[NA/L] for any integer L, we 
observe, therefore, that if k is replaced by k—1 and N by NA, while leaving 
NV unchanged*, the results (15) to (20) will still hold for the present case. 


Thus we obtain 


k-1 
sup|S(A) — NV|< J] [M,(R; —1) + (2R;—1)]< 
t=1 (25) 





k-1 
R,;-—2 
< (InN) J] (3% if N>every R;, 
( ) IT( In R; ) ery 
and, similarly, 
k-1 k-1 (R; —1)? 
J<4- OTT (M+ 3) (R12 <4 YT] SR, nN + 31nRi)*, (26) 


whence, if N > max (R}), for example, 


T<(2- —(k— -0 FT ms ee) an N24 1) (27) 


We see that these results are the same as (5) and (6), with 


(R; — 1)* 
eer “Ea og (28) 


a= TT nk) 


t=1 





5. Some further remarks 


The result (5) matches VAN DER CoRPUT’s inequality (4) when k=2. It 
confirms HAMMERSLEY’S conjecture. The result (6) is supplementary to this, 
and of comparable strength. There is still, unfortunately, a considerable gap 
between (6) and Rotn’s lower bound (3). There are evidently three possible 
causes for this. (i) My results may not be as strong as possible. (ii) There may 
be a different sequence from that of VAN DER CorRPUT and HAMMERSLEY, for 
which a stronger result than mine may be deduced (just as the sequence (10) 
gives a stronger result, namely (6), than the result (24) for the sequence I propose 
in (11)). (iii) RoTH’s result may be too weak. Of these, I have considered only 
the first two possibilities. 





* The V in (17) to (20) must be interpreted as A, A, ... A,_, and multiplied by 
NA; but we now wish to consider V = A A, A,... Ax_,; thus NV must be left un- 
changed. 
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Consider (i). The result (18) is exact. The step to (19) involves an increase 
in value by a factor of the order of /7;(1+M;-"), which would not appreciably 
close the gap between (3) and (6). From (20) to (24), an approximation of the 
same order is again all that occurs. It would seem, therefore, that the step from 
(19) to (20), im which the sign and magnitude of y= (8 — {N/IT,R?"}) is concerned, 
is the only one where an improvement might be hoped for, to lead to a strengthen- 
ing of the result (6). 

In the absence of any detailed information about the behavior of y, we may 
observe that, if y is distributed symmetrically about zero (that is, if +yp are 
equally frequent), then the sums on the right of (19) will average out to a constant 
multiple of [J7;(M;+2)]—# times the product on the right of (20). Thus sup 
|S(A)—NV| and |/J will each be of the order of (InN)*” for the sequence (41), 
and so of the order of (InN)“*~?” for the Hammersley sequence. We may thus 
tentatively conjecture that positive constants B,, C), exist, such that 


sup|S(A) — NV| <C,(InNO-”) (29) 
and 
J < B,(nN)*}, (30) 


for the Hammersley sequence, with similar results for the sequence (11). 


Consider now (ii). I have considered the effect, for given N and a one- 
dimensional sequence, yz(1), Pe(2), ---, Pr(N), of type (11), of weighting points 
according to the number of digits in the corresponding integer. The (R” — R”—') 
points generated by m-digit integers are given a weight W,. Then, if we put 
W=1, Wy+1=0, and S,,=W, — W,,., (m=0, 1, 2,..., M), we get the relations 





M M 
Sm = Wo— Wy ii=1, S,R™=N, 
m=0 m=0 
M (34) 
|S (A) —NV| =| >, S,,{A R™} — | 
m=0 
These lead to 
M M m-1 
J=3 2 Smt £2 2D Su Sat +$R™). (32) 


The unweighted case gives S,,=0 (m=O, 1, 2,..., M—1), Sy=1, whence 


N=R™, J=3-. (33) 
Note that, owing to the selection of N as a power of R, the InN factor has been 
eliminated. Thus we see that (33) 1s an tmprovement on VAN DER CORPUT’S r¢sti/i 
for N=R™ only. It can be shown that the corresponding k-dimensional resu): 
simply multiplies the (InN)?* factor by (k —1)*/k*. (Only in the case of k=1 
does this essentially alter the nature of J.) If we minimise (32), subject to the 
first equation in (31), we obtain the result that J ~ 7, so that unequal weighting 
does not seem to be worth-while. 


6. Computation of sequences 


The following scheme is suitable for computing the successive members of 
the sequence, 


Pr(1), Pr(2), Pr(3), +--+» PR(N),---, (34) 
7* 
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using only the operations available in digital computers. The operations listed 
below are carried out consecutively, unless a “jump” instruction intervenes. 


Put p=0 
(I) Put x=1-— 9 GQivsis”  dnssdadabvatducaseed sb iseuntecebhedeeeane Enter 
Put y=1/R 


Put new y = old y/R 
Jump to (LD) one eeseeeeeeee 


(III) Put new g = old 9+ (R+1)y— 





The three-step loop starting at (II) finds the most significant (t.e. leftmost) 
digit, in the radix-R expression for the current fraction gy, which is less than 
(R—1). The next value of y is then computed by (III). For each successive 
value of y, the routine is entered at (I) and left at (III). 





I wish to thank Dr. J. M. HAMMERSLEY, for bringing his paper [1] to my notice, 
and the United Kingdom Atomic Energy Authority, for a research grant. 
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A Numerical Method for the Solution of a Parabolic System 
By 


Jim Douctas jr. 


1. Introduction 


The object of this paper is to introduce a finite difference method for the 
approximate solution of a parabolic system of differential equations that arises 
in the description of two-phase flow of incompressible fluids in a multi-dimensional 
porous medium. The differential system is 


V-(aVu)+V-(BVv) =0, 


1.1 
oe V-(BVu)+V-(aVv) =p 22, 

where the coefficients «, 8, and y are functions of x, y, and v. The coefficients 
satisfy the relations 


MZ=«a(x,y,v) =m>0, 
(1.2) | B(x, y, v)/a(x, y,v)| SB<1, 
y(x,y,v) =m>0. 


In addition, they are expected to be twice continuously differentiable. In the 
analysis to follow, the coefficients could also depend explicitly on ¢, or on more 
space variables, but not on wu. Physically, « represents the sum of the pressures 
in the two phases, v the difference between the pressures in the two phases, 
« the sum of the mobilities of the two phases, B the difference of these mobilities, 
and y the rate of change of the saturation of one phase with respect to v. See 
[8] for a derivation of (1.1) and applications of the method to be given below 
to specific physical problems, including comparison with laboratory data. 

In order to avoid extremely small time steps, implicit difference methods 
are preferred. Because of the nature of the boundary conditions that arise in 
certain problems, alternating direction methods for the transient system analogous 
to those having proved useful for one parabolic equation [5,9] do not appear 
universally applicable. As a result, a backward difference system, 


A(a,, A U,, +1) +A(B,4 Vi.41) =0, 


Va+1 —v, 


“ A(By A Uy41) + Alp AV 41) = t=, 








This research was supported in part by the Air Force Office of Scientific Research 
under Contract Number AF 49(638)632. 











92 Jim Dovuctas jr.: 


is suggested. The symbols in (1.3) can be defined as below: 
x;=tAx, y,=j4y, t,=nAt, 
Vi, = Vin = V(X, Hj bn), 
O, = 2(%;,9;,Vizn),  ebe., 
(1.4) 41/2 =% (% 41/2, ¥j>%(Visaj,nt+ Visa), 
A, (a, 4, U,,41) = [os 4s2(Uiss ints = U; j,n+1) = 
— %12(U; 5 nt — U;-1,j,n41) ]/(4%)?, 
A(a,, AU, +41) =A, (a, A,U,.41) + A, (a, 4, U,41) - 





For convenience, assume 4x= Ay. 

The convergence of the solution of the difference system (1.3) to that of 
the differential system (1.1) will be proved for a slight simplification of (1.1) 
in Section 2. In Section 4 an alternating direction iteration method for solving 
the system of linear equations that occur at each time step will be introduced, 
and under strong restrictions the convergence of this iteration procedure will be 
established. 

2. Convergence 


Ht is a commonly accepted point of view that lower order terms in a differential 
operator do not affect the convergence or divergence of a difference analogue 
of the differential operator. Of course, such a principle as this is purely heuristic ; 
however, let us apply the principle to (1.1). Then (1.1) becomes 


aM4u+fpAv=0, 


(24) pAut+adv=y%, 


where A is the Laplace operator. Our convergence proof will be given for (2.1), 
which may be simplified to 


Au+adAv=0, 
(2.2) aAu+Av=b<, \a(x,y,v)|SB<14, 
b(x, y,v) =m>0, 
by division by «. The difference system reduces to 
(42 + A}) Una +44 (43 + 45) Vass =0, 


—V"n 
a, (A2 + A3) Uy 4, + (42 + 43) Vay =, HR , 


(2.3) 


where 4? indicates the centered, divided second difference. 
Simple eliminations give 


dv 


(1—a*)Av=b—, 


(2.4) 


Vat — Va 
At x 


(1 — a3) (4? + Aj) Via =b, 
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So far, we have not specified anything about either the region in which 
(2.2) is to be satisfied or the boundary conditions to be assigned to u and v. 
Assume, for the purposes of this analysis, that the region is a finite cylinder in 
(x, y, )-space with the base a polygon of segments of the lines x= x; and y= yj 
for all choices of 4x considered. Also, assume that values of u and v are assigned 
on the boundary and v initially, and assume that u and v are four times boundedly 
differentiable in the region. 


Under the above assumptions, it is easy to see by the argument of [4] that 
(2.5) max|v—V|=O((4x)?+At), OS#tST, 


for any T as Ax and 4? tend to zero independently. Thus, V is a satisfactory 
approximation to v uniformly in the region. 


The analysis of the error in u is more difficult. First, let us derive a difference 
equation for the error. Let 


6, =u, — U,, 
(2.6) | 


&, =v, — V,. 


Under the assumptions made above, 


(A? + A?) Uns +4(%, Y,Un41) (A? + Aj) Vn41 =O ((4x)?) . 
Thus, 


(2.7) (43 +A}) Un 41 +4(%, ¥, Op) (AE + A5) 0441 =O ((4x)* +42). 
Subtracting the first equation of (2.3) from (2.7) and simplifying, we obtain 
(42 + 45) Butt + On (42 +4) Ents + Cy &, =O ((4x)? + At), c, =O(1). 

As ¢,=0((4x)?+ 42), 
(2.8) (43 + A3) 5,41 +44 (4? + A3) e444 =O ((4x)? + Ad). 
Recall that 

A? (f:8:) =f, 458; + Asha 4e8i41 + Af 428i + 8 4Eh, 


A,f;=(f;—fi-1)/4*. 


where 


Thus, 


2 2 Bes 2 aa 2 
Oy, M5 Eng = Ae (An Ent) — Ay G41), 4s €541,j,n+1 A,4;;,4, Ej j,n+1 £n41424y- 





Now, 
a(x;, 94 Viin) — @ (44-1, 94 Vi-w in) 
A,4;;,= (¥i, Hj, Viin oe 1» Vj» “i-1,5,0 
0a da 
= Ge tay 4sMim 
7] da 
= 4 7 7 ~y (4,5, — 4&5 n]- 
Assume that 
At 


(2.9) =constant>0, p22. 


(4%) p 
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Then, 4,¢,=O(4x), and each of the terms 4,a,4,¢, is O(4x). Similarly, 
Aia,=O(1) and ¢,42a,=O((4x)*). Consequently, 

(2.10) (43 A}) (6,41 + Oy En+1) = 0 (4x) 


at interior grid points and 6,,,—¢,,,;=0 at boundary points. Then, by an a 
priori bound for solutions of elliptic difference equations due to Bers [J, 2], 


(2.14) Ont1 + 4n ngs =O(42). 
As En41 =O ((4x)?*), 
(2.42) 8,4, =0(Az). 


Thus, under the assumption outlined, convergence of the solution of the difference 
system (2.3) to the solution of the differential system (2.2) is assured. 


3. A Generalization 
The analysis given in the last section can be extended in a straightforward 
manner to the system 


3.1) | meee 


cM4u+Av+d=ey,, 


where a, c, and e depend on v and the independent variables and b and d depend 
on v, v,, vy, and the independent variables. It is necessary to assume that 


jac|\ SB<i, 
(3.2) | 5. +|5,,| +|4,,| +|d,,| SM<oo, 
e=m>o0, 


along with the usual smoothness conditions. In the difference equation, v, and 
v, are replaced by symmetric differences at the new time level in evaluating 
b and d. 


4. An Alternating Direction Iteration Procedure 


The system (1.3) is clearly implicit; fortunately, the algebraic equations at 
each time step are linear, since «, 8, y are evaluated with known data. The 


algebraic problem is of the generic form 
A(a AU) + A(BAV) =0, 
se A(BAU) + A(a AV) - V =f, |B|S Ba, y>0o. 


For the argument below, it is not necessary that B<1, although we treat it as 
if it were. Consider the following iteration procedure: 


A,(aA,U**** + BA, V+) + A, (ad, U*+f.4,V" 
= A,(a UFt+2 +B yeti? _ 4 Ut — B V*), 
A,(BA,U**** +04, V*+*) + A,(B4,U* +04,V*) — 2 yeti 
=f +A,(BU**** + @ V+ — BUt—a V4), 


(4.2a) 
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and 


A, (aA, U*+¥? + BA, V*t¥) + A, (aA, t+? + BA, V+) 
=A, (a Ut! 4 BVA+t — @ UF+U2 _ Byhtiey | 

A,(BA,U*+*? +A, V***) + A, (B.A, U** +04, Vit) — 2 yet 
=f +A, (B U*t! + Vit? — BUFtN — 4 YR) | 


(4.2b) | 





This iteration method has been used at each time step, using the values 
of U and V at the old time level as initial guesses, for several problems. It proved 
to be about ten times as efficient as over-relaxation for a region of approximately 
100 grid points. The over-relaxation was not done with an optimum choice of 
the over-relaxation parameter, as the coefficients changed at each time step. 

It is well-known [3, 10] that alternating direction iteration methods are hard 
to analyze. We shall have to restrict «, 8, and y to be constants and the region 
to be a rectangle, which for convenience will be assumed to be the unit square. 
We may absorb the constant y into the time ¢ by a linear change of variables; 
consequently, let y=1. Then, the change of dependent variables 


ahi 


(4.3) Breet 


simplifies (4.2) to 


AEB" + AFB = AG? — BP), 
SA | A ABT at BPM — 0M) — AO +1, 
and 

A? peta +A? pit} — en sae prtnay | 
(4.4b) AR g*t* + AS gt 4+ —____ a= ile By = (Bp t* — ag*t!) — A, (q**+? —q**") +f. 
Let 

p= —P, 

‘ (ee-3 
(4.5) Fae 


Then, the convergence of the iteration method is equivalent to p*-0, g*—>0. 
Obviously, ~* and g* satisfy the homogeneous form of (4.4) with vanishing 
boundary values. The eigenfunctions of this form of (4.4) are easily seen to be 
of the form 


(4.6) 
Let 


{ p* =P,sinnmxsinany, 


gq’ =Q,sinamxsinany. 


a 
© = Fie 7 
A 
(4.7) X= Tas sin® i e, 





4 2 mnAx 


= (Az)? sin 3 
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By straightforward substitution, we see that 











(4.8 (A, +X) Py4ij2=(A, — Y) P,, 
8a 
(A, +X + C) Qn +1/2 — £ C Pys12 1 (A, uaa Y) Or, 
and 
(4.80) (A, + Y) Py = (Ay — X) Phiye, 
; (A,+Y+C) Orta — 2 C Phas = (An —X) Qr+1/2- 
Thus, 
P,,. = (Ar—X) (4r—Y) 
k+1 (Ay +X) (Ant Y) k? 
De (Ax—X) (Axn— Y) 
— | P= Ca +X+C) (At ¥ +E) Pht 


-+- BC (Ay,—X) (Ap — Y) 

a | (Ay+X) (Ap +Y+C) (A, +X4+C) 
(A,—X) (Ax—Y) , 

(Ap +X) (Ap + Y) (A, +Y+C) | * 


It can be shown that, if A,=A for all k, P* and Q* tend to zero for all m, 
n=1,...,(4x)+—1. Thus, the choice of a constant iteration parameter leads 
to a convergent process. However, this is by no means the most efficient choice. 
Note that P, is independent of Q,; moreover, it has exactly the same recursion 
relation as the recursion relation arising when Laplace’s equation is treated by 
an alternating direction method [9]. Thus, the repeated use of a certain finite 
geometric sequence for {A,} will assure that the total number of iterations necess- 
ary to reduce all ((4x)"+— 1)? of the P, to a preassigned level will be no moxie 
than O(—log4x), as opposed to O((4x)-*) for the choice of constant A,. See 
[6,9]. Let us see that Q, can also be made small in but a slightly larger number 
of iterations. 

Clearly, for the ratio of P,,, to P, to be less than one in magnitude, A, must 
lie between the smallest and largest values of X. Thus, 


4 
(4%) 


p 





of 








9g 2Ax 
2 


(4.10) min A, =min X = ~ 





; sin 
as Ax—>0, and 


BC/a BC/a B 
(4.14) afte < A <|£| <1. 








Then, if 
_ | (Ar—X) (Ar— Y) | 
aie * | (Ap+X) (Ag+ Y) |’ 
(4.13) | Onss| Sel Ox| +204] F- 
Now, 
k-1 
Consequently 


k 
(4.15) | Qn+1| S oe| Qe! + 2| Po| IT 0;. 
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It is easy to see that the solution of the recursion relation 


hy arbitrary, 








(4.16) 
hy s1 = Onhy + cil Qi» 
is 7" 
k—1 
(4.17) hy = (hg + Ch) TT Q;. 
i=0 
Thus, 
k-1 
(4.18) | Qu] S (1 Qo] + 2| Pol 4) IT a: 
By [6, 9], we may choose {A,} such that 
Anis=Ay, 
r logaAx/2 
(4.19) , tog ((1—84)/(1 +48) 
s-—1 
0,<6<1, all mand xn. 
i=0 
Then, 
k-1 [=| Lay 
IT 0;<4 SEM. tg 
i=0 
and 
k 
eae | 
(4.20) 1Qze| S (| Qo] +2] Plz) ds , all mand n. 


In the above analysis, P, and Q, depend on mand 0; 1.¢., P, =P, mn» On = Qe, min: 


Let us use 
(4x)-1 


(4.24) p ; (PE mn + QO m»)| =E, 


as the norm of the error. Then, 


k 
Eis » oo + (| Qo, m,n| > 2k|Po m,nl)*} sls _ 


(4.22) < sls -) > {(4 + 8k?) ae +2 Q5, m,n 
< gargtls =) Ei. 
Thus, 
. & 
(4.23) E,S3kd* Ey. 


Let the iteration be carried out until the ratio of E, to E, is reduced to a given 
level: 


ae 
Zk* 65 Se, 


h*—1 


3(k*—1)6 * >. 


(4.24) 


Let z=k*/s and solve 


—, 
z 6 ae 
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for z. It is known [10] that the best choice of 6 is y2-1 for the alternating 
direction process for the Laplace difference equation. Since one of our relations 
is the Laplace difference equation, let us also use this value. Then, 


z— 1.14logz=1— 1.14 log ‘ 
As logz/z<e™, 1.14logz< 2/2 and 
2<2— 2.28 log ¢/3s 


=O (log - + log log 3) ‘ 
Thus, 
(4.25) k* =O ((tog a + log log =z) log ss) ; 


Since the difference system (4.2) is a pair of 7-diagonal linear equations, Gaussian 
elimination applied to this system requires O ((4x)~*) calculations per iteration. 
Thus, for the constant coefficient case in a rectangle, no more than 


O ((4 x)~* log ae (log “ -+ log log =z)} 


calculations are required to obtain a satisfactory solution of (4.1). 


In applications to the physical examples of [8], it was observed that only 
one cycle of s iterations was needed for most time steps of the transient problem. 
Except for the first two or three time steps, no more than two cycles were ever 
required. 
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On the Numerical Determination 
of the Best Approximations in the Chebyshev Sense 
By 
L. VEIDINGER 


The problem of the Chebyshev approximation, in case of real functions of 
a real variable, can be generally formulated as follows (see [1], p. 1): let there 
be given the set of points H of the real line, and the real functions f(x) and 
F(x; A,, Ag, ..-,A,) of the real variable x, which are bounded on this set. We 
have to determine the values Af, A, ..., AX of the parameters /,, A,,..., A, for 
which the maximum deviation between F(x; A,, Ag, ..., A,) and f(x) 


E (Ay, Ag, -++1 An) = max | f(x) — F(x; Ay, Ag, ..-, Ay)| 


is the least, provided such parameter values exist and can be given uniquely. 
Henceforward we shall restrict ourselves to the case when the set of points H 
is a finite and closed interval [a, b], the function /(x) is continuous throughout 
this interval, and 


F(x; A, Ae, sila A.) = A, 9 (x) + Ay Po (*) om + Ay Pn (*) (1) 


where the real functions 9,(x), y2(x), ..., ~,(%) are linearly independent and 
continuous on [a, b] and form a CHEBYSHEV’s system of functions over this 
interval. This means that any function of the form (1) (hence these functions 
will be called polynomials), which does not vanish identically in [a,b] may 
not have more than m—1 zeros in this interval. (It is agreed that inside the 
interval [a,b] each zero of the polynomial F(x; /,, A,,...,A4,) at which the 
polynomial does not change sign, must be considered as double.) As it is 
well-known, (see for example [1], p. 74), the polynomial of best approximation 
F(x) =F (x; Af, ag, ..., AS) has the property that the difference 


A(x) = f(x) — F(x) 
takes the value 
: E = max |A(x)| 


at least at +1 successive points of the interval [a,b] with alternating signs. 
This property characterizes the polynomial F(x)=F(x; Af, Az, ..., Af) uniquely. 

The polynomial F(x) can be determined exactly only in the simplest cases, 
therefore we have to apply numerical methods for approximation. The problem 
to determine numerically the best approximations in the Chebyshev sense became 
of great practical importance as soon as the modern digital computers appeared. 
If we want to compute the values of a transcendental function by means of a 
computer it is most convenient to use a rational function approximating the 








100 L. VEIDINGER: 


given function over a sufficiently large interval rather than interpolating in the 
table of the function-values. The best approximations in the Chebyshev sense 
are probably the most suitable ones for this purpose. It must be noticed that 
it is not always practical to choose the form (1) for these best approximations 
(see [2], p.95—114), but the theory of the Chebyshev approximation in the 
case of more general functions F(x; A,, A,,...,4,) is still in the rough. (For 
results in this direction see [3] and [4].) 

One of the simplest numerical methods to determine the polynomial of best 
approximation in the Chebyshev sense is the so-called exchange method (see, 
for example, [10], p. 224). The exchange method is as follows. 


Let us choose the points 
Ky < HES <M 


in [a, 6] as we please. It is easy to determine the polynomial F, (x)= Ajg, (x) + 
A32(x)+---+Alq,(x) which approximates the function (x) best in the Cheby- 
shev sense on the set of the points x (7=1, 2,...,m+1), if we.start from the 
fact that the difference 

A, (x) = f(x) — F(x) 


at the points x}, attains equal absolute values with alternating signs*: 
A, (%}) = 1(9)) — BA P(A) = LE (— A, G=12--.m+1) (2) 


where E,=0 (see for example [7], p. 176). The coefficients Aj, 4}, ..., Al can be 
evaluated by solving the former system of linear equations. (The determinant 
of the coefficients is other than zero, as we shall see below.) 

Let x! be a point where |4,(x)| takes its maximum in [a, b]. If the point x1 
coincides with one of the points x}, then F, (x)= F(x): otherwise we may replace, 
as it is easy to see, one of the points x} by the point x! so that at the points 


KE< MEM <M, 
of the so obtained new point set the function 4,(x) takes its values with alter- 
nating signs. Let us construct the polynomial F,(x)=Aj@, (x) + A3q_(*)+---+ 
42 p,(x) approximating the function f(x) best on the points x4 (j=1, 2, ...,m+4); 
then 
A, (xj) = f(%;) —F(x})=+8,(—1) ((=1,2,...,."+1) 
where E,=0, and the method is to be continued in the same way. 


The idea of the convergence proof of REMEs, given in [5], p. 338—339, may 
be applied to the exchange method, and in the particular case when 9, (x)=1, 
P(x) =, ..., P,(x)=%"—*, gives that for all sufficiently large values of v: 


E—E,Sq(E — E,_,) 


where gq is independent of vy and 0<q<1 (see [5], p. 339); that is, the difference 
E—E, tends to zero “at least as rapidly” as g’. The proof of REMEs can be 
easily extended to the general case too. 





* Here and henceforward we may ascribe an arbitrary sign to the zero function- 
values. 
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A more elaborate algorithm for the determination of the polynomials of best 
approximation was discovered by REMEs (see [5]; this is the so-called second 
algorithm of REMES) and was described in various modified forms by many 
authors in the late years (see [2], [3], [6], [8], and [9]). The principle of the 
Remes method is similar to the principle of the exchange algorithm. In order 
to emphasize the analogy, we shall not alter the notations introduced above. 


Let the points 
n< nh< or eas 


be arbitrarily chosen in [a, 6], and let the points 
MS Mesos < 24, (v= 2, },...) 


be formed by the »—1-th step. Let us construct the polynomials F,(x) resp. 
F,(x) to these points in the same way as in the enemy poseens, but pass from the 
points x7~* to the points x7 in gt way. Let z{~*, 23°", ..., 2-1 be namely zeros 
of the function A,_,(*)= fix) — _, (x) for which wise tc x51 (¢=4, 2,..., #) 
holds; it is obvious that such zeros a tered Let us elect from each of the sub-intervals 
[a, Ze 251], (k=1,2,..., —1) and [z?-*, 6] a point, where A,_,(x) takes 
its maximum or minimum in the orresponding sub-interval, depending on the 
positive or negative sign of the function-values 4,_,(x;~*), 4,_,(%%53) resp. 

A,_1(%,31). We denote by y;~* the point elected from the j-th sub-interval, 
going from left to right. If the function 4,_,(x) takes its absolute maximum 
in [a, 6] at one of the points y;~* then let 43=y;~* (j=1,2 ...,m+4) other- 
wise we replace as in the exchange method, one of the points y;~* by such a 
point x’~* where |4,_,(x)| takes its maximum in [a, b] and now we denote the 
j-th point of the so obtained new point set by x;; going from left to right. 


We shall show in the remaining part of this paper that under certain conditions 
in the Remes algorithm: 
E — E,=O(E — E, _,)*. (3) 


That is, in this case the difference E — E, tends to zero “at least as rapidly” 
as r?”, where 0<r<1 and 7 is independent of ». “ 


The fact that E —E,->0 follows from the general proof of NOvoDVORSKII 
and PINSKER (see [7]). We can deduce immediately from their proof that the 
polynomials F(x) converge uniformly to the polynomial F(x) in (a, 6]: moreover, 
denoting by x; (j=1, 2, ..., #+ 4) some of the accumulation points of the sequence 
x; (v=4, 2, ...) then 

A(x;) = + E(—1)! 


and (in case of 7<m+1) %;<%4, 

Henceforward we shall suppose the functions /’(x) and gj; (x) (¢=4, 2, ..., ») 
to be continuous in the closed interval [a, 6] and continuously derivable (once) 
in the open interval (a, b). We also suppose the function |4(x)| to take its maxi- 
mum in [a,b] exactly at +41 points. (Every maximum point falling inside 
[a,b], where A”’(x) vanishes and every maximum point coinciding with either 
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of the end points of the interval where 4’(x) vanishes must be considered as 
double!) * 

It follows immediately from our restriction for the function 4(x) that the 
sequences of the points x; (j= 1, 2, ..., m+ 4) are convergent. For the polynomials 
F,(x) tend uniformly to the polynomial F(x) in [a,b], Af—>A¥ (i=, 2,..., n) 
and hence it follows that the functions F’(x) and 4’ (x) converge-uniformly to the 
functions F’(x) resp. A’(x) in [a, b]; furthermore if a<a’<b’<b the functions 
F-"(x) and A?’(x) also tend uniformly to the functions F’’(x) resp. A’’(x) in [a’, 6’). 

To prove the inequality (3) we need the following lemma: 

If 

Uy < Ug <*< Mya 


are arbitrary points of [a, b], the determinants 


V1 (4) Pa(%) «+. Py (4) 


v, (4j~1) Po(%j~1) --- Pn (M1) 


(j =1,2,...," +14) 
YP (Mj 41) Po(4j41) --- Pn (Mj 41) 








Pi (Un+1) Pa(Mn+a) «++ Pn (Mn41) 
are other than zero and of equal sign. 


In order to prove the lemma, let us consider for 1=1, 2,..., the polynomial 


1 (4) a (%) --- Py (4) 


v1 (4:1) Po(%;-1) --- Pn (M1) 
®; ;1(u) =| 9, (4) Po(™) —... Py (4) 
1 (U;42) Po(%j+2) --- Pn(%+2) 








Pi(Un+1) Po(Mn41) «++ Pn (Mn+a) 
= my? py (m) + Mgt t* palm) +--+ + Rt” Pn (u)- 


The polynomial ®; ;,,;(“) may have at most —1 zeros in [a,b], because not 
each of the coefficients ui**?, uy**?, ..., ub**?, is equal to zero (see [1], p. 69). 
The polynomial vanishes at the points 4,,..., #1, 449, +--+» 4% 41, and so the 
values ®; ;,,(u;)=D;,, and ®; ;,,(4;,,;)=D,; are not zero and of equal sign. 
(From the lemma it follows immediately that, the determinant of the system 
of equations (2), and similarly, for y=2,3,..., the determinant of the system 


He) — LHe) =£E(— 1) G=A,2, +4) 


do not vanish.) 





* We were unable to give a general criterion to state the above mentioned prop- 
erty of the function 4(*) from the given function f(*) and the ‘“‘basic polynomials’”’ 
yp; (*) (j= 7 ae 93-04 n). 
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From the construction of the points x; it follows that if the point x, is an 
interior point of the interval [a, b], then 


A, -1(%) = (49) — xa 9i(%})=0. (j= 41,2,...," +14) (4) 


holds. If the point x; coincides with either a or b, and » is sufficiently large, 
then %3= a". For, if j in the cases x,=a or x,,,=6, we could find a number 9, 
such that x zt (resp. %n+1) would be interior points of the interval [a, b] for all 
values of ash then, since xj->a; (x,,,,->6) and the functions A’ (x) tend uni- 
formly to the function 4’(x) in [a, b], it would follow that A’(a) =0 (resp. A’(b) =0). 
But this is in contradiction to our condition for the function A(x). From the 
above discussion it follows that 


A,-1(%) = (a) — Dae p(x) 
= +E,_,(—1)'— Ses g —x-)s  (j=1,2,...,8 +1) (5) 


where, if x= x;~*, & is between the points x; and x}, and for all sufficiently 
large values of Y, sg A" (&)= =+(—1)'*?; if f= 27-1 , then & is an arbitrary point 
in (a,b). If E(x)= ¥ w(x) + Ber (s)-+-- 42 9a) is the best approximation 
of the function 4,_,(%) on the point set {xj}, and A, (x) =A,(x) — E(x), then 


A, 4) = A,1(4) — Dhl) =£E (1 GHA). 6 
In fact, if 


A, (2) = (a) — [B-1(4) + Ea] = 4 E(— 1) G=4,2,....8+4) 


where E’ > E,, then hence and from the relations 
A, (x) =H (9) —E()=4E,(—1))  G§=1,2,....0+1) 


it would follow that the polynomial F,_,(x)+ F(x) —E(x), not identically zero 
in [a, b], would have at least m zeros in this interval. 
The relations (6) yield a system of %+1 linear equations for the unknowns 
1,Ag,---,A, and E,. By the lemma, the determinant of the coefficients in the 
system of equations is other than zero. Therefore, we may express the unknown 
E, from the system of linear equations (6) by means of CRAMER’S rule: 


Pos |A3| |4,—1 (44)|+ [45] 14,—1 (43)|+ °** + 14n4s| 14, i(¥n+a)l (7) 
4 |A¥|+ |A3|+ °° + |Ahaal 





where 
V1 (x1) Pa(%i) «+» Pn (*1) 


(4a) ala) «++ Pal4}a) | 
tile) Pa (%j41) +++ Pn (%}41) 








nl n+1) P2(%,41) -- + Pn (%41) 
Numer. Math. Bd. 2 8 
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Substituting (5) into (7), and using the fact that the determinants A; and their 
limits the determinants 


1 (%;) Po(%) — -++ Pu(%) 


4 (%;-1) P2(%j—-1) «+» Pn (%j—1) 


1 (%;41) Po(%j41) «+» Pn(% 41) 





P1(%n41) Po(%n41) «+» Pn (Xn41) 





in virtue of the lemma, differ from zero it follows that 
at} ” 1 att ” 
6, 2, |4,-1(&)| (4, — 4; PSE,—E£,15 22, 14,1 (&)| (x5 — a3~")*, (8) 
7= I= 


where the constants c, and c, are independent of ». Naturally the numerical 
values of the constants depend on m, on the function f(x) and on the ‘“‘basic 
polynomials” ;(x) (¢=1,2,...,#). The same holds for the later constants 
C3, +++, Cz too. As 


4,()) =H) — La pile) =+E(— 1, G=1,2,040+1) 
hence and from (5) we get 


14” v— y j 4, i » I— y 
(i — 84) oa) = 4 (EE, (— 1 + SP a age ig 


7 
(j=1,2,...,8+4). 


ims 


t 


The relations (9) yield a system of +1 linear equations for the unknowns 
x; —H-* (i=1, 2,...,). We express the unknown 4;—{~* from the first m 
equations of the system (9) by means of CRAMER’S rule, and get 


Re — a a Mt Mi a Fine gaan) (10) 








An+1 
where M%,, Mis, ..., Mj, are the minors of the elements of the i-th column in 
the determinant Aj,,, respectively. From (8) and (9) it follows that 
y Gi ; 
| Q5| <(1 + x} Z,- E,_1) = ¢,(E, — E,_;) j= 1, 2, ee A+ 1). 


Since the minors M{,, Mjz,..., Mj, are less in absolute value than a positive 
bound independent of v, (as the “basic polynomials” ;(x) are bounded) and 
the determinants A’,,, and A,,, differ from zero, we obtain from (10) that 


|A;—A*| Sc (E,—E,_,) (¢=1,2,...,m). (11) 


The functions |qg;(x)| are all less than a common bound (independent of »), 
therefore (4) and (11) yield 


| 4; (%5)| =|4,(43) — 4,_1(%))| = \> (4; — AF") gi (%})| Se5(E, — E,-1). 


t=1 


1 
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From (4), and by means of the mean value theorem, we get the inequalities 
A’ (x7+2) — At (at)| =|” (0%) (a8? — 22)| < 06 (E, — E, 1) (12) 


where the point ¢7 is between x; and xj** in the case when aj** +x, and Cs 
is an arbitrary point in (a, 6) if atten gf If the limit point x; of the Sequence 
x; will lie inside the interval [a, b], then, from our restriction for A(x), A’’(x,) +0, 
but the sequence of the functions Aj’(x) (y=1, 2,...) converges uniformly to 
A’’(x) in a sufficiently small neighbourhood of x,;, and so from (12) it follows, 
for all sufficiently large values of » that 


|x;*" ans x5 | S c,(E, — E,-_,)- 


If x,=a or x, ,,=5, then, as we have seen before, if v is large enough, x}= x3+!=a 
resp. *,,1=,14=b. We deduce from the above statements and from the 
inequality 


1 


n+1 ati 
o4 2 1Ay (B5*)| (9f** — 29)* S Ey 44 — Ey Sea 214, ()| (24? — 25), 
a 2 


analogous with (8), the relation 
E,41 a E, = C9 (E, aa E,_,)* 
and hence we immediately get the inequality (3). 
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Introduction 


Background. After the publication!.? of a preliminary report on the algorith- 
mic language ALGOL, as prepared at a conference in Ziirich in 1958, much 
interest in the ALGOL language developed. 


As a result of an informal meeting held at Mainz in November 1958, about 
forty interested persons from several European countries held an ALGOL im- 
plementation conference in Copenhagen in February 1959. A ‘“‘hardware group” 
was formed for working cooperatively right down to the level of the paper tape 
code. This conference also led to the publication by Regnecentralen, Copenhagen, 
of an ALGOL Bulletin, edited by PETER NAuR, which served as a forum for 
further discussion. During the June 1959 ICIP Conference in Paris several 
meetings, both formal and informal ones, were held. These meetings revealed 
some misunderstandings as to the intent of the group which was primarily re- 
sponsible for the formulation of the language, but at the same time made it 
clear that there exists a wide appreciation of the effort involved. As a result of 
the discussions it was decided to hold an international meeting in January 1960 
for improving the ALGOL language and preparing a final report. At a European 
ALGOL Conference in Paris in November 1959 which was attended by about 
fifty people, seven European representatives were selected to attend the January 
1960 Conference, and they represent the following organisations: Association 
Francaise de Calcul, British Computer Society, Gesellschaft fiir Angewandte 
Mathematik und Mechanik, and Nederlands Rekenmachine Genootschap. The 
seven representatives held a final preparatory meeting at Mainz in December 
1959. 

Meanwhile, in the United States, anyone who wished to suggest changes or 
corrections to ALGOL was requested to send his comments to the ACM Com- 
munications where they were published. These comments then became the basis 
of consideration for changes in the ALGOL language. Both the SHARE and USE 
organisations established ALGOL working groups and both organisations were 





1 Preliminary report — International Algebraic Language, Comm. Assoc. Comp. 
Mach. 1, No. 12 (1958), 8. 

2 Report on the Algorithmic Language ALGOL by the ACM Committee on Pro- 
gramming Languages and the GAMM Committee on Programming, edited by A. J. 
PerR.tis and K. SAMELSoN, Numerische Mathematik Bd. 1, S. 41—60 (1959). 
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represented on the ACM Committee on Programming Languages. The ACM 
Committee met in Washington in November 1959 and considered all comments 
on ALGOL that had been sent to the ACM Communications. Also, seven repre- 
sentatives were selected to attend the January 1960 international conference. 
These seven representatives held a final preparatory meeting in Boston in De- 
cember 1959. P 


January 1960 Conference. The thirteen representatives!, from Denmark, 
England, France, Germany, Holland, Switzerland, and the United States, con- 
ferred in Paris from January 11 to 16, 1960. 

Prior to this meeting a completely new draft report was worked out from the 
preliminary report and the recommendations of the preparatory meetings by 
PETER NAUvR and the Conference adopted this new form as the basis for its report. 
The Conference then proceeded to work for agreement on each item of the report. 
The present report represents the union of the committee’s concepts and the 
intersection of its agreements. 

As with the preliminary ALGOL report, three different levels of language 
are recognized, namely a Reference Language, a Publication Language and 
several Hardware Representations. 


Reference Language 


1. It is the working language of the committee. 

2. It is the defining language. 

3. The characters are determined by ease of mutual understanding and not 
by any computer limitations, coders notation, or pure mathematical notation. 

4. It is the basic reference and guide for compiler builders. 

5. It is the guide for all hardware representations. 

6. It is the guide for transliterating from publication language to any locally 
appropriate hardware representations. 

7. The main publications of the ALGOL language itself will use the reference 
representation. 


Publication Language 


1. The publication language admits variations of the reference language ac- 
cording to usage of printing and handwriting (e.g., subscripts, spaces, exponents, 
Greek letters). . 

2. It is used for stating and communicating processes. 

3. The characters to be used may be different in different countries, but 
univocal correspondence with reference representation must be secured. 


Hardware Representations 


1. Each one of these is a condensation of the reference language enforced by 
the limited number of characters on standard input equipment. 





1 WILLIAM TURANSKI of the American group was killed by an automobile just 
prior to the January 1960 Conference. 
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2. Each one of these uses the character set of a particular computer and is the 
language accepted by a translator for that computer. 

3. Each one of these must be accompanied by a special set of rules for trans- 
literating from publication or reference language. 


For transliteration between the reference language and a language suitable 
for publications, among others, the following rules are recommended. 


Reference Language Publication Language 

Subscript brackets [ ] Lowering of the line between the brackets and 
removal of the brackets. 

Exponentiation + Raising of the exponent. 

Parentheses ( ) Any form of parentheses, brackets, braces. 

Basis of ten 49 Raising of the ten and of the following integral 
number, inserting of the intended multiplica- 
tion sign. 


Description of the reference language 


Was sich tiberhaupt sagen 148t, 146t sich 
klar sagen; und wovon man nicht reden 
kann, dartiber mu68 man schweigen. 


Lupwic WITTGENSTEIN 
1. Structure of the language 


As stated in the introduction, the algorithmic language has three different 
kinds of representations—reference, hardware, and publication—and the develop- 
ment described in the sequel is in terms of the reference representation. This 
means that all objects defined within the language are represented by a given set 
of symbols—and it is only in the choice of symbols that the other two represen- 
tations may differ. Structure and content must be the same for all representations. 


The purpose of the algorithmic language is to describe computational processes. 
The basic concept used for the description of calculating rules is the well known 
arithmetic expression containing as constituents numbers, variables, and func- 
tions. From such expressions are compounded, by applying rules of arithmetic 
composition, self-contained units of the language—explicit formulae—called 
assignment statements. 


To show the flow of computational processes, certain nonarithmetic state- 
ments and statement clauses are added which may describe e.g., alternatives, or 
iterative repetitions of computing statements. Since it is necessary for the func- 
tion of these statements that one statement refers to another, statements may be 
provided with labels. Sequences of statements may be combined into compound 
statements by insertion of statement brackets. 


Statements are supported by declarations which are not themselves comput- 
ing instructions, but inform the translator of the existence and certain properties 
of objects appearing in statements, such as the class of numbers taken on as 
values by a variable, the dimension of an array of numbers or even the set of 
rules defining a function. Each declaration is attached to and valid for one 
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compound statement. A compound statement which includes declarations is 
called a block. 


A program is a self-contained compound statement, i.e. a compound state- 
ment which is not contained within another compound statement and which 
makes no use of other compound statements not contained within it. 


In the sequel the syntax and semantics of the language will be given?. 


1.1. Formalism for syntactic description. 


The syntax will be described with the aid of metalinguistic formulae*. Their 
interpretation is best explained by an example: 


<ab> ::= (| [| <ab>( | <ab) <a 


Sequences of characters enclosed in the bracket <)> represent metalinguistic 
variables whose values are sequences of symbols. The marks ::= and | (the 
latter with the meaning of or) are metalinguistic connectives. Any mark in a 
formula, which is not a variable or a connective, denotes itself (or the class of 
marks which are similar to it). Juxtaposition of marks and/or variables in a 
formula signifies juxtaposition of the sequences denoted. Thus the formula above 
gives a recursive rule for the formation of values of the variable <ab). It indi- 
cates that <ab)> may have the value ( or | or that given some legitimate value of 
<ab)>, another may be formed by following it with the character ( or by follow- 
ing it with some value of the variable <d>. If the values of <d)> are the decimal 
digits, some values of <ab) are: 

[(((2(37( 

(12345( 


((( 
[86 


In order to facilitate the study the symbols used for distinguishing the metalin- 
guistic variables (i.e. the sequences of characters appearing within the brackets 
<)> as ab in the above example) have been chosen to be words describing approxi- 
mately the nature of the corresponding variable. Where words which have 
appeared in this manner are used elsewhere in the text they will refer to the 
corresponding syntactic definition. In addition some formulae have been given 
in more than one place. 


Definition: 


<empty>:: = 
(i.e. the null string of symbols). 





1 Whenever the precision of arithmetic is stated as being in general not specified, 
or the outcome of a certain process is said to be undefined, this is to be interpreted 
in the sense that a program only fully defines a computational process if the accom- 
panying information specifies the precision assumed, the kind of arithmetic assumed, 
and the course of action to be taken in all such cases as may occur during the exe- 
cution of the computation. 

2 Cf. J. W. Backus, The syntax and semantics of the proposed international 
algebraic language of the Ziirich ACM-GAMM conference. ICIP Paris, June 1959. 
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2. Basic symbols, identifiers, numbers, and strings. 
Basic concepts 


The reference language is built up from the following basic symbols: 
<basic symbol): : = <letter)| <digit >| (logical value)| (delimiter) 


2.1. Letters 
<letter> :: = alb|c|d|e|f|g|h|1|j||2|m|n|0|p|q|r|s|¢||v| | x] y|2| 
A|B|C|D|E|F|G|A|Z|J|K|L|M|N|O| P| Q|R|S|7|U|V|W |X| ¥|Z 
This alphabet may arbitrarily be restricted, or extended with any other distinc- 
tive character (i.e. character not coinciding with any digit, logical value or 
delimiter). 
Letters do not have individual meaning. They are used for forming identifiers 
and strings! (cf. sections 2.4. Identifiers, 2.6. Strings). 


2.2.1. Digits 
<digit> :: = 0| 1| 2| 3| 4| 5|6|7| 8|9 
Digits are used for forming numbers, identifiers, and strings. 


2.2.2. Logical values 
<logical value :: = true| false 
The logical values have a fixed obvious meaning. 


2.3. Delimiters 
<delimiter) :: = <operator)| <separator)| <bracket)| <declarator)| <specificator) 
<operator) :: = arithmetic operator)| «relational operator)| <logical operator) | 
<sequential operator) 


<arithmetic operator) ::= + | —|x|/|+|t 

<relational operator) ::= <|<| =|2|>|+ 

<logical operator) :: = =|>|V|A|— 

<sequential operator) :: = go to| if| then| else | for| do? 
<separator>::=,|.|19|:|;|:= | u |step| until] while | comment 
<bracket> ::= (|) | []]] ‘|’ | begin] end 


<declarator) :: = own | Boolean | integer | real | array | switch | procedure 
<specificator) :: = string | label | value 

Delimiters have a fixed meaning which for the most part is obvious, or else 
will be given at the appropriate place in the sequel. 

Typographical features such as blank space or change to a new line have 
no significance in the reference language. They may, however, be used freely 
for facilitating reading. 





1 It should be particularly noted that throughout the reference language under- 
lining (for typographical reasons synonymously bold type. Publishers remark) is 
used for defining independent basic symbols (see sections 2.2.2 and 2.3). These 
are understood to have no relation to the individual letters of which they are 
composed. 

2 do is used in for statements. It has no relation whatsoever to the do of the 
preliminary report, which is not included in ALGOL 60. 
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For the purpose of including text among the symbols of a program the fol- 
lowing “‘comment”’ conventions hold: 


The sequence of basic symbols: is equivalent with 
; comment <any sequence not containing ;); . 
begin comment <any sequence not containing ;); begin 
end <any sequence not containing end or ; or else) end 


By equivalence is here meant that any of the three symbols shown in the right 
hand column may, in any occurrence outside of strings, be replaced by any se- 
quence of symbols of the structure shown in the same line of the left hand column 
without any effect on the action of the program. 


2.4. Identifiers 


2.4.1. Syntax. 
<identifier) :: = <letter)| <identifier) <letter>| <identifier) (digit > 
2.4.2. Examples q 
Soup 
Vi17a 
a34kTMNs 
MARILYN 


2.4.3. Semantics. Identifiers have no inherent meaning, but serve for the iden- 
tification of simple variables, arrays, labels, switches, and procedures. They 
may be chosen freely (cf. however section 3.2.4. Standard functions). 

The same identifier cannot be used to denote two different quantities except 
when these quantities have disjoint scopes as defined by the declarations of the 
program (cf. section 2.7. Quantities, kinds and scopes and section 5. Declara- 
tions). 


2.5. Numbers 
2.5.1. Syntax. 
<unsigned integer) :: = (digit>| (unsigned integer) <digit> 
<integer> :: = (unsigned integer) | -+ «unsigned integer)>| — «unsigned integer) 
<decimal fraction) :: =. unsigned integer) 
<exponent part) :: = 49 <integer> 
<decimal number) :: = <unsigned integer)>| (decimal fraction)| 
<unsigned integer) «decimal fraction) 
<unsigned number) :: = (decimal number)| exponent part)| 
_ <decimal number) <exponent part)» 
<number) :: = (unsigned number)|-} (unsigned number)| — «unsigned number) 


2.5.2. Examples. 0 — 200.084 —.083)—02 
177 407.4318 rel 
5384 9.34,)+10 all 
+0.7300 Q—4 +5 


2.5.3. Semantics. Decimal numbers have their conventional meaning. The ex- 
ponent part is a scale factor expressed as an integral power of 10. 

2.5.4. Types. Integers are of type integer. All other numbers are of type real 
(cf. section 5.1. Type declarations). 
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2.6. Strings 

2.6.1. Syntax. 

<proper string» :: = <any sequence of basic symbols not containing ‘ or ’| 
<empty> 


<open string) :: =<proper string>|‘<open string>’| <open string) <open string) 
<string> :: =‘<open string»’ 
2.6.2. Examples. 

‘d5k,, agi ‘TIA =/ ? Tt’ 

‘.. This uts ua u‘string’’ 


2.6.3. Semantics. In order to enable the language to handle arbitrary sequences 
of basic symbols the string quotes ‘ and ’ are introduced. The symbol u denotes 
a space. It has no significance outside strings. 

Strings are used as actual parameters of procedures (cf. sections 3.2. Function 
designators and 4.7. Procedure statements). 


2.7. Quantities, kinds and scopes 

The following kinds of quantities are distinguished: simple variables, arrays, 
labels, switches, and procedures. 

The scope of a quantity is the set of statements in which the declaration for 
the identifier associated with that quantity is valid, or, for labels, the set of 
statements which may have the statement in which the label occurs as their 
successor. 


2.8. Values and types 


A value is an ordered set of numbers (special case: a single number), an 
ordered set of logical values (special case: a single logical value), or a label. 


Certain of the syntactic units are said to possess values. These values will in 
general change during the execution of the program. The values of expressions 
and their constituents are defined in section 3. The value of an array identifier 
is the ordered set of values of the corresponding array of subscripted variables 
(cf. section 3.1.4.1). , 

The various ‘‘types” (integer, real, Boolean) basically denote properties of 
values. The types associated with syntactic units refer to the values of these units. 


3. Expressions 


In the language the primary constituents of the programs describing algorith- 
mic processes are arithmetic, Boolean, and designational, expressions. Consti- 
tuents of these expressions, except for certain delimiters, are logical values, 
numbers, variables, function designators, and elementary arithmetic, relational, 
logical, and sequential, operators. Since the syntactic definition of both variables 
and function designators contains expressions, the definition of expressions, and 
their constituents, is necessarily recursive. 


<expression) :: = (arithmetic expression) | «Boolean expression) | 
<designational expression) 
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3.1. Variables 


3.1.1. Syntax. 

«variable identifier) :: = <identifier) 

<simple variable» :: = <variable identifier) 
<subscript expression) :: = <arithmetic expression) 


<subscript list : : = (subscript expression) |<subscript list), (subscript expression) 
<array identifier) :: = <identifier) 

<subscripted variable) :: = array identifier) [<subscript list] 

<variable> :: = ¢simple variable>|<subscripted variable) 


3.1.2. Examples. epsilon 
det A 
al7 
Q[7, 2] 
x[sin(n x pi/2), Q[3, n, 4] 


3.1.3. Semantics. A variable is a designation given to a single value. This value 
may be used in expressions for forming other values and may be changed at will 
by means of assignment statements (section 4.2). The type of the value of a 
particular variable is defined in the declaration for the variable itself (cf. section 
5.1. Type declarations) or for the corresponding array identifier (cf. section 5.2. 
Array declarations). 


3.1.4. Subscripts. 3.1.4.1. Subscripted variables designate values which are 
components of multidimensional arrays (cf. section 5.2. Array declarations). Each 
arithmetic expression of the subscript list occupies one subscript position of the 
subscripted variable, and is called a subscript. The complete list of subscripts 
is enclosed in the subscript brackets []. The array component referred to by a 
subscripted variable is specified by the actual numerical value of its subscripts 
(cf. section 3.3. Arithmetic expressions). 


3.1.4.2. Each subscript position acts like a variable of type integer and the 
evaluation of the subscript is understood to be equivalent to an assignment to 
this fictitious variable (cf. section 4.2.4). The value of the subscripted variable 
is defined only if the value of the subscript expression is within the subscript 
bounds of the array (cf. section 5.2. Array declarations). 


3.2. Function designators 


3.2.1. Syntax. 
<procedure identifier) ::= <identifier) 
<actual parameter) ::= <string> | <expression) | <array identifier) | 


<switch identifier) | <procedure identifier) 
<letter string) ::= <letter) | <letter string) <letter) 
<parameter delimiter) ::=, | ) <letter string) :( 
<actual parameter list) ::= actual parameter) | 

<actual parameter list) <parameter delimiter) <actual parameter) 
<actual parameter part) ::= <empty) | (actual parameter list) 
<function designator) ::= <procedure identifier) <actual parameter part) 
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3.2.2. Examples. sin (a — b) 
J(v+s,n) 
R 


S(s — 5) Temperature : (T) Pressure : (P) 

Compile (‘: =’) Stack : (Q) 
3.2.3. Semantics. Function designators define single numerical or logical values, 
which result through the application of given sets of rules defined by a procedure 
declaration (cf. section 5.4. Procedure declarations) to fixed sets of actual para- 
meters. The rules governing specification of actual parameters are given in 
section 4.7. Procedure statements. Not every procedure declaration defines the 
value of a function designator. 
3.2.4. Standard functions. Certain identifiers should be reserved for the standard 
functions of analysis, which will be expressed as procedures. It is recommended 
that this reserved list should contain: 
abs (E) for the modulus (absolute value) of the value of the expression E 
sign (E) for the sign of the value of E(+ 1 for E>0, 0 for E=0, — 1 for E< 0) 
sqrt (E) for the square root of the value of E 
sin(E) for the sine of the value of E 
cos (E) for the cosine of the value of E 
arctan (E) for the principal value of the arctangent of the value of E 
In (E) for the natural logarithm of the value of E 
exp(E) for the exponential function of the value of E (e*) 
These functions are all understood to operate indifferently on arguments both 
of type real and integer. They will all yield values of type real, except for sign (E) 
which will have values of type integer. In a particular representation these 
functions may be available without explicit declarations (cf. section 5. Declara- 
tions). 
3.2.5. Transfer functions. It is understood that transfer functions between any 
pair of quantities and expressions may be defined. Among the standard functions 
it is recommended that there be one, namely 


entier (E), 
which “‘transfers” an expression of real type to one of integer type, and assigns 
to it the value which is the largest integer not greater than the value of E. 


3.3. Arithmetic expressions 

3.3.4. Syntax. adding operator) ::=+ | — 

<multiplying operator) ::= x | / | + 

<primary> ::= <unsigned number) | variable) | <function designator) | 
(<arithmetic expression) 

<factor> ::= <primary) | <factor> + <primary) 

<term) ::= <factor> | <term)> <multiplying operator) <factor) 

<simple arithmetic expression) ::= <term) | <adding operator) <term) | 
<simple arithmetic expression) <adding operator) <term) 

<if clause) ::= if <Boolean expression) then 

<arithmetic expression) ::= <simple arithmetic expression) | 
<if clause> <simple arithmetic expression) else <arithmetic expression) 
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3.3.2. Examples. 
Primaries: 

7.394,,—8 

sum 

w[t+ 2, 8] 

cos (y + 2x3) 

(a — 3/y+vut 8) 


Factors: 
omega 
sum * cos(y + z x3) 
7.394,9—8t w[i+2,8]t(a—3/y+vuut 8) 


Terms: 


U 
omega X sum t cos (y + z X3)/7.394,) — 8+ w[i + 2, 8] t (a — 3/y +-v ut 8) 


Simple arithmetic expression: 
U — Y «+ omega x sum + cos (y+2z X3)/7.394,.—8 t w[i + 2, 8] t (a — 3/y+vut 8) 


Arithmetic expressions: 


wxu—Q(S+Cu)t2 

if g>0 then S+3xQ/A else 2xS+3xq 

ifa<0 then U-+ V else if axb> 17 then U/V else if k + y then V/U else 0 
a X sin (omega xt) 

0.57912 xa[N x(N — 1)/2, 0] 

(A xarctan(y) +Z) + (7 + Q) 

if g then  — 1 else n 

if a< 0 then A/B else if b= 0 then B/A else z 


3.3.3. Semantics. An arithmetic expression is a rule for computing a numerical 
value. In case of simple arithmetic expressions this value is obtained by exe- 
cuting the indicated arithmetic operations on the actual numerical values of the 
primaries of the expression, as explained in detail in section 3.3.4. below. The 
actual numerical value of a primary is obvious in the case of numbers. For 
variables it is the current value (assigned last in the dynamic sense), and for 
function designators it is the value arising from the computing rules defining 
the procedure (cf. section 5.4. Procedure declarations) when applied to the current 
values of the procedure parameters given in the expression. Finally, for arithmetic 
expressions enclosed in parentheses the value must through a recursive analysis 
be expressed in terms of the values of primaries of the other three kinds. 


In the more general arithmetic expressions, which include if clauses, one out 
of several simple arithmetic expressions is selected on the basis of the actual 
values of the Boolean expressions (cf. section 3.4. Boolean expressions). This 
selection is made as follows: The Boolean expressions of the if clauses are eva- 
luated one by one in sequence from left to right until one having the value true 
is found. The value of the arithmetic expression is then the value of the first 
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arithmetic expression following this Boolean (the largest arithmetic expression 
found in this position is understood). The construction: 


else <simple arithmetic expression) 


is equivalent to the construction: 
else if true then <simple arithmetic expression) 





ne rears net ete mca 


3.3.4. Operators and types. Apart from the Boolean expressions of if clauses, 
the constituents of simple arithmetic expressions must be of types real or integer | 
(cf. section 5.1. Type declarations). The meaning of the basic operators and the i 
types of the expressions to which they lead are given by the following rules: 


3.3.4.4. The operators +, —, and x have the conventional meaning (addition, 
subtraction, and multiplication). The type of the expression will be integer if 
both of the operands are of integer type, otherwise real. 


3.3.4.2. The operations <term)/<factor>) and <term)— <factor>) both denote 
division, to be understood as a multiplication of the term by the reciprocal of 
the factor with due regard to the rules of precedence (cf. section 3.3.5). Thus for 
example 

alb x7/(p — q) xvjs 
means 

((((@ x (6)) x7) x((6 — g)*)) xv) x(s™) 
The operator / is defined for all four combinations of types real and integer 
and will yield results of real type in any case. The operator ~ is defined only 
for two operands both of type integer and will yield a result of type integer 
defined as follows: 

a — b = sign(a/b) x entter (abs (a/b)) 
(cf. sections 3.2.4 and 3.2.5). 


3.3.4.3. The operation <factor> + <primary) denotes exponentiation, where the 
factor is the base and the primary is the exponent. Thus for example 


2tntk means (2")* 
while 
2t(ntm) means 2” 


Writing i for a number of integer type, 7 for a number of real type, and a for ) 
a number of either integer or real type, the result is given by the following rules: | 


ati Ift>0: axax---xXa (t times), of the same type as a. | 
Ifi=0, if a+0: 1, of the same type as a. 
if a=0O: undefined. 
Ifi<0, ifa+0: I1/(axax--- Xa) (the denominator has 
+ factors), of type real. | 
if a=0: undefined. ! 
atr If a>0: exp(rx/n(a)), of type real. 
Ifa=0, if r>0: 0.0, of type real. 
if y= 0: undefined. 
If a<0: always undefined. 








a ee ik re RA OR i 
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3.3.5. Precedence of operators. The sequence of operations within one expression 
is generally from left to right, with the following additional rules: 

3.3.5.1. According to the syntax given in section 3.3.1 the following rules of 
precedence hold: 


first: x 
second: x/— 
third: +— 


3.3.5.2. The expression between a left parenthesis and the matching right paren- 
thesis is evaluated by itself and this value is used in subsequent calculations. 
Consequently the desired order of execution of operations within an expression 
can always be arranged by appropriate positioning of parentheses. 


3.3.6. Arsthmetics of real quantities. Numbers and variables of type real must 
be interpreted in the sense of numerical analysis, i.e. as entities defined inherently 
with only a finite accuracy. Similarly, the possibility of the occurrence of a 
finite deviation from the mathematically defined result in any arithmetic ex- 
pression is explicitly understood. No exact arithmetic will be specified, however, 
and it is indeed understood that different hardware representations may evaluate 
arithmetic expressions differently. The control of the possible consequences of 
such differences must be carried out by the methods of numerical analysis. 
This control must be considered a part of the process to be described, and will 
therefore be expressed in terms of the language itself. 


3.4. Boolean expressions 


3.4.1. Syntax. 
<relational operator) ::= <| <| =|2=|>|+ 
<relation) ::= <arithmetic expression) <relational operator) 


<arithmetic expression 
<Boolean primary) ::= <logical value) | <variable> | <function designator) | 

<relation) | («Boolean expression) 
<Boolean secondary) ::= <Boolean primary) | — <Boolean primary) 
<Boolean factor) ::= <Boolean secondary) { 

<Boolean factor) A <Boolean secondary» 
<Boolean term) ::= <Boolean factor) | <Boolean term) V <Boolean factor) 
<implication) ::= Boolean term) | <implication> > <Boolean term) 
<simple Boolean) ::= <implication> | <simple Boolean) = <implication) 
<Boolean expression) ::= <simple Boolean) | 

_ <if clause) <simple Boolean) else <Boolean expression) 


3.4.2. Examples. x = — 2 

Y>VVz<q 

a+b>—d5Az-d>qt2 

pAQGV x+y 

g=—aAbA-—cVdVer-f 

if k<1thens>w elsehsc 

if if if a then b else c then d else / then g else h< hk 
3.4.3. Semantics. A Boolean expression is a rule for computing a logical value. 
The principles of evaluation are entirely analogous to those given for arithmetic 
expressions in section 3.3.3. 





ak i iii 
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3.4.4. Types. Variables and function designators entered as Boolean primaries 
must be declared Boolean (cf. section 5.1. Type declarations and section 5.4.4. 
Values of function designators). 


3.4.5. The operators. Relations take on the value true whenever the corresponding 
relation is satisfied for the expressions involved, otherwise false. 

The meaning of the logical operators — (not), A (and), V (or), > (implies), 
and = (equivalent), is given by the following function table. 





b41 false false true true 
b2 false true false true 
— 51 true true false false 
b1Ab2 false false false true 
b1V 62 false true true true 
b1362 true true false true 
bi =bd2 true false false true 














3.4.6. Precedence of operators. The sequence of operations within one expression 
is generally from left to right, with the following additional rules: 


3.4.6.1. According to the syntax given in section 3.4.1. the following rules of 
precedence hold: 

first: arithmetic expressions according to section 3.3.5. 

second: << S$ = 2> + 

third: — 

fourth: A 

fifth: V 

sixth: > 

seventh: = 


3.4.6.2. The use of parentheses will be interpreted in the sense given in section 
3.3.5.2. 


3.5. Designational expressions 


3.5.4. Syntax. 

<label> ::= <identifier> | «unsigned integer) 

<switch identifier) ::= <identifier) 

<switch designator) ::= <switch identifier) [<subscript expression >] 

<simple designational expression) ::= <label) | <switch designator) | 
(<designational expression) 

<designational expression) ::= <simple designational expression) | 
<if clause) <simple designational expression) else <designational expression > 


3.5.2 Examples. 17 


p9 

Choose [n — 1] 

Town [if y < 0 then N else N + 1) 

if Ab<c then 17 else g[if w< 0 then 2 else n] 
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3.5.3. Semantics. A designational expression is a rule for obtaining a label of 
a statement (cf. section 4. Statements). Again the principle of the evaluation 
is entirely analogous to that of arithmetic expressions (section 3.3.3). In the 
general case the Boolean expressions of the if clauses will select a simple designa- 
tional expression. If this is a label the desired result is already found. A switch 
designator refers to the corresponding switch declaration (cf. section 5.3. Switch 
declarations) and by the actual numerical value of its subscript expression 
selects one of the designational expressions listed in the switch declaration by 
counting these from left to right. Since the designational expression thus 
selected may again be a switch designator this evaluation is obviously a re- 
cursive process. 


3.5.4. The subscript expression. The evaluation of the subscript expression is 
analogous to that of subscripted variables (cf. section 3.1.4.2). The value of a 
switch designator is defined only if the subscript expression assumes one of the 
positive values J, 2, 3, ..., m, where ” is the number of entries in the switch list. 


3.5.5. Unsigned integers as labels. Unsigned integers used as labels have the 
property that leading zeroes do not affect their meaning, e.g. 00217 denotes the 
same label as 217. 


4. Statements 


The units of operation within the language are called statements. They will 
normally be executed consecutively as written. However, this sequence of 
operations may be broken by go to statements, which define their successor 
explicitly, and shortened by conditional statements, which may cause certain 
statements to be skipped. 


In order to make it possible to define a specific dynamic succession, statements 
may be provided with labels. 


Since sequences of statements may be grouped together into compound 
statements and blocks the definition of statement must necessarily be recursive. 
Also since declarations, described in section 5, enter fundamentally into the 
syntactic structure, the syntactic definition of statements must suppose de- 
clarations to be already defined. 


4.1. Compound statements and blocks 
4.1.1. Syntax. 
<unlabelled basic statement) ::= <assignment statement) | <go to statement) | 
<dummy statement) | <procedure statement) 
<basic statement) ::= <unlabelled basic statement) | <label) : <basic statement) 
<unconditional statement) ::= <basic statement) |<for statement) | 
<compound statement) | <block> 
<statement)> ::= <unconditional statement) | <conditional statement) 
<compound tail) ::= <statement) end | <statement) ; <compound tail) 
<block head) ::= begin <declaration> | <block head) ; <declaration) 
<unlabelled compound) ::= begin <compound tail) 
<unlabelled block) ::= <block head) ; <compound tail) 


Numer. Math, Bd. 2 9 
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<compound statement) ::= <unlabelled compound) | <label) : «compound 
statement» 
<block> ::= <unlabelled block) | <label> : <block) 
This syntax may be illustrated as follows: Denoting arbitrary statements, de- 
clarations, and labels, by the letters S, D, and L, respectively, the basic syntactic 
units take the forms: 
Compound statement: 
L: L:... beginS ;S;...S; Send 
Block: 
L: L: ... begin D; D;..D;S;S;...S; Send 
It should be kept in mind that each of the statements S may again be a complete 
compound statement or block. 


4.1.2. Examples. Basic statements: 
ai=p+q 
go to Naples 
Start: Continue: W := 7.993 


Compound statement: 
begin x:=0; for y:=1 step J until n do x:=x+ Aly]; 
if x >q then go to STOP else if x > w — 2 then go to S; 
Aw: St: W:=x-+ bob end 


Block: 
Q: begin integer 7, k; real w; 
for 1:= 1 step J until m do 
for k:=i-+ 1 step J until m do 
begin w:= A[1, k]; 
Aft, k]:= A[k, 7]; 
A[k,t]:=w end fori and k 
end block Q 


4.1.3. Semantics. Every block automatically introduces a new level of nomen- 
clature. This is realized as follows: Any identifier occurring within the block 
may through a suitable declaration (cf. section 5. Declarations) be specified to 
be local to the block in question. This means (a) that the entity represented 
by this identifier inside the block has no existence outside it and (b) that any 
entity represented by this identifier outside the block is completely inaccessible 
inside the block. 


Identifiers (except those representing labels) occurring within a block and 
not being declared to this block will be non-local to it, i.e. will represent the same 
entity inside the block and in the level immediately outside it. The exception 
to this rule is presented by labels, which are local to the block in which they 
occur. 


Since a statement of a block may again itself be a block the concepts local 
and non-local to a block must be understood recursively. Thus an identifier, 
which is non-local to a block A, may or may not be non-local to the block B in 
which A is one statement. 
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4.2. Assignment statements 
4.2.1. Syntax. 
<left part) ::= <variable> := 
<left part list) ::= <left part) | <left part list) <left part) 
<assignment statement) ::= <left part list) <arithmetic expression) | 
<left part list) <Boolean expression) 
4.2.2. Examples. s:= [0]: =n:=n+1+s 
n:=n+1 
A:= B/C —v—qxS 
s[v,k + 2] := 3 — arctan(s x zeta) 
V:=Q>YAZ 
4.2.3. Semantics. Assignment statements serve for assigning the value of an 
expression to one or several variables. The process will in the general case be 
understood to take place in three steps as follows: 
4.2.3.1. Any subscript expressions occurring in the left part variables are eva- 
luated in sequence from left to right. 


4.2.3.2. The expression of the statement is evaluated. 


4.2.3.3. The value of the expression is assigned to all the left part variables, 
with any subscript expressions having values as evaluated in step 4.2.3.1. 
4.2.4. Types. All variables of a left part list must be of the same declared type. 
If the variables are Boolean the expression must likewise be Boolean. If the 
variables are of type real or integer the expression must be arithmetic. If the 
type of the arithmetic expression differs from that of the variables, appropriate 
transfer functions are understood to be automatically invoked. For transfer 
from real to integer type the transfer function is understood to yield a result 
equivalent to 

entier (E + 0.5) 
where E is the value of the expression. 


4.3. Go to statements 
4.3.1. Syntax. 
<go to statement) ::= go to <designational expression 
4.3.2. Examples. 

go to 8 

go to exit[n+ 1] 

go to Town [if y < 0 then N else N + 1) 

go to if Ab<c then /7 else g[if w < 0 then 2 else n} 
4.3.3. Semantics. A go to statement interrupts the normal sequence of operations, 
defined by the write-up of statements, by defining its successor explicitly by the 
value of a designational expression. Thus the next statement to be executed 
will be the one having this value as its label. 
4.3.4. Restriction. Since labels are inherently local, no go to statements can lead 
from outside into a block. 
4.3.5. Go to an undefined switch designator. A go to statement is equivalent 
to a dummy statement if the designational expression is a switch designator 
whose value is undefined. 


9* 
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4.4. Dummy statements 
4.4.1. Syntax. 
<dummy statement) ::= <empty) 
4.4.2. Examples. 


} &. 
begin ....; John: end 


4.4.3. Semantics. A dummy statement executes no operation. It may serve to 
place a label. 


4.5. Conditional statements 


4.5.1. Syntax. 

<if clause> :: = if <Boolean expression) then 

<unconditional statement) :: = <basic statement) | <for statement) | 
<compound statement)| <block> 

<if statement) :: = <if clause) (unconditional statement | 
<label>: <if statement > 

<conditional statement) :: = <if statement>| <if statement) else (statement) 


4.5.2. Examples. if x>0 then n:=n+1 
if v>u then V:¢:=n-+m else goto R 
if s<0V PSQ then AA: begin if g<v then a: =v/s 
else y: =2 x a end else if v>s then a: =v—gq 
else if vu>s— 1then go to S 


4.5.3. Semantics. Conditional statements cause certain statements to be executed 
or skipped depending on the running values of specified Boolean expressions. 


4.5.3.1. If statement. The unconditional statement of an if statement will be 
executed if the Boolean expression of the if clause is true. Otherwise it will be 
skipped and the operation will be continued with the next statement. 


4.5.3.2. Conditional statement. According to the syntax two different forms of 
conditional statements are possible. These may be illustrated as follows: 


if B1 then S1 else if B2 then S2 else S3; S4 
and 
if B1 then S1/ else if B2 then S2 else if B3 then S3; S4 


Here B1 to B3 are Boolean expressions, while S1 to S3 are unconditional state- 
ments. S4 is the statement following the complete conditional statement. 


The execution of a conditional statement may be described as follows: The 
Boolean expressions of the if clauses are evaluated one after the other in sequence 
from left to right until one yielding the value true is found. Then the uncondi- 
tional statement following this Boolean is executed. Unless this statement de- 
fines its successor explicitly the next statement to be executed will be S4, i.e. 
the statement following the complete conditional statement. Thus the effect 
of the delimiter else may be described by saying that it defines the successor of 
the statement it follows to be the statement following the complete conditional 


statement. 
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The construction 

else «unconditional statement) 
is equivalent to 

else if true then <unconditional statement) 

If none of the Boolean expressions of the if clauses is true, the effect of the 
whole conditional statement will be equivalent to that of a dummy statement. 


For further explanation the following picture may be useful: 


t t \ 
if B1 then S1 else if B2 then S2 else S3; S4 
+S. co aa t 


7 Bi false. B2 false 


4.5.4. Go to into a conditional statement. The effect of a go to statement leading 
into a conditional statement follows directly from the above explanation of the 
effect of else. 


4.6. For statements 
4.6.1. Syntax. 
<for list element >:: = <arithmetic expression) | 
<arithmetic expression) step <arithmetic expression) until 
<arithmetic expression) | 
<arithmetic expression) while «Boolean expression) 
<for list>:: = <for list element)| <for list), <for list element) 
<for clause>:: = for <variable>: = <for list) do 
<for statement) :: = <for clause) (statement) | 
<label>: <for statement) 
4.6.2. Examples. for g: =I step s until n do A[q]: = B[q] 
for k:=1, V1 x 2 while V1<N do 
for j3:=1+G, L, 1 step J until N, C+D do A[k,j|: =B[k,7] 
4.6.3. Semantics. A for clause causes the statement S which it precedes to be re- 
peatedly executed zero or more times. In addition it performs a sequence of 
assignments to its controlled variable. The process may be visualized by means 
of the following picture: 


} 
Initialize; test; statement S; advance; successor 
M bats 


for list exhausted 


In this picture the word initialize means: perform the first assignment of the for 
clause. Advance means: perform the next assignment of the for clause. Test 
determines if the last assignment has been done. If so the execution continues 
with the successor of the for statement. If not the statement following the for 
clause is executed. 

4.6.4. The for list elements. The for list gives a rule for obtaining the values 
which are consecutively assigned to the controlled variable. This sequence of 
values is obtained from the for list elements by taking these one by one in the 
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order in which they are written. The sequence of values generated by each of 
the three species of for list elements and the corresponding execution of the state- 
ment S are given by the following rules: 

4.6.4.1. Arithmetic expression. This element gives rise to one value, namely 
the value of the given arithmetic expression as calculated immediately before 
the corresponding execution of the statement S. 


4.6.4.2. Step-until-element. An element of the form A step B until C, where A, 
B, and C, are arithmetic expressions, gives rise to an execution which may be 
described most concisely in terms of additional ALGOL statements as follows: 
V:=A; 
L1: if (V —C) x sign(B)>0 then go to Element exhausted ; 
Statement S; 
V:=V+B; 
go to LI; 
where V is the controlled variable of the for clause and Element exhausted 
points to the evaluation according to the next element in the for list, or if the 
step-until-element is the last of the list, to the next statement in the program. 
4.6.4.3. While-element. The execution governed by a for list element of the 
form E while F, where E is an arithmetic and F a Boolean expression, is most 
concisely described in terms of additional ALGOL statements as follows: 
L3:V:=E; 
if — F then go to Element exhausted; 
Statement S; 
go to L3; 
where the notation is the same as in 4.6.4.2. above. 
4.6.5. The value of the controlled variable upon exit. Upon exit out of the state- 
ment S (supposed to be compound) through a go to statement the value of the 
controlled variable will be the same as it was immediately preceding the execution 
of the go to statement. 
If the exit is due to exhaustion of the for list, on the other hand, the value 
of the controlled variable is undefined after the exit. 
4.6.6. Go to leading into a for statement. 
The effect of a go to statement, outside a for statement, which refers to a 
label within the for statement, is undefined. 


4.7. Procedure statements 


4.7.1. Syntax. 
<actual parameter): : = (string>| <expression)| (array identifier) | 
<switch identifier) | «procedure identifier) 
<letter string>:: = <letter)| <letter string) <letter) 
<parameter delimeter):: =, |) (letter string): ( 
<actual parameter list>:: = (actual parameter) | 
<actual parameter list) <parameter delimeter) «actual parameter) 
<actual parameter part):: = <empty>| (actual parameter list)) 
<procedure statement): : = <procedure identifier) <actual parameter part)» 
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4.7.2. Examples. Spur (A) Order: (7) Result to: (V) 

Transpose (W, v + 1) 

Absmax (A, N, M, Yy, I, K) 

Innerproduct (A(t, P, u], B[P], 10, P, Y) 
These examples correspond to examples given in section 5.4.2. 


4.7.3. Semantics. A procedure statement serves to invoke (call for) the execution 
of a procedure body (cf. section 5.4. procedure declarations). Where the proce- 
dure body is a statement written in ALGOL the effect of this execution will be 
equivalent to the effect of performing the following operations on the program: 


4.7.3.1. Value assignment (call by value). All formal parameters quoted in the 
value part of the procedure declaration heading are assigned the values (cf. 
section 2.8. Values and types) of the corresponding actual parameters, these 
assignments being considered as being performed explicitly before entering the 
procedure body. These formal parameters will subsequently be treated as local 
to the procedure body. 


4.7.3.2. Name replacement (call by name). Any formal parameter not quoted 
in the value list is replaced, throughout the procedure body, by the corresponding 
actual parameter, after enclosing this latter in parentheses wherever syntactically 
possible. Possible conflicts between identifiers inserted through this process and 
other identifiers already present within the procedure body will be avoided by 
suitable systematic changes of the formal or local identifiers involved. 


4.7.3.3. Body replacement and execution. Finally the procedure body, modified 
as above, is inserted in place of the procedure statement and executed. 


4.7.4. Actual-formal correspondence. The correspondence between the actual 
parameters of the procedure statement and the formal parameters of the proce- 
dure heading is established as follows: The actual parameter list of the procedure 
statement must have the same number of entries as the formal parameter list 
of the procedure declaration heading. The correspondence is obtained by taking 
the entries of these two lists in the same order. 


4.7.5. Restrictions. For a procedure statement to be defined it is evidently 
necessary that the operations on the procedure body defined in sections 4.7.3.1. 
and 4.7.3.2 lead to a correct ALGOL statement. 

This poses the restriction on any procedure statement that the kind and type 
of each actual parameter be compatible with the kind and type of the correspond- 
ing formal parameter. Some important particular cases of this general rule are 
the following: 


4.7.5.1. Strings cannot occur as actual parameters in procedure statements 
calling procedure declarations having ALGOL 60 statements as their bodies 
(cf. section 4.7.8). 


4.7.5.2. A formal parameter which occurs as a left part variable in an assignment 
statement within the procedure body and which is not called by value can only 
correspond to an actual parameter which is a variable (special case of expression). 
4.7.5.3. A formal parameter which is used within the procedure body as an 
array identifier can only correspond to an actual parameter which is an array 
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identifier of an array of the same dimensions. In addition if the formal para- 
meter is called by value the local array created during the call will have the same 
subscript bounds as the actual array. 

4.7.5.4. A formal parameter which is called by value cannot in general correspond 
to a switch identifier or a procedure identifier, because these latter do not possess 
values (the exception is the procedure identifier of a procedure declaration which 
has an empty formal parameter part (cf. section 5.4.1) and which defines the 
value of a function designator (cf. section 5.4.4). This procedure identifier is in 
itself a complete expression). 

4.7.5.5. Any formal parameter may have restrictions on the type of the corre- 
sponding actual parameter associated with it (these restrictions may, or may not, 
be given through specifications in the procedure heading). In the procedure 
statement such restrictions must evidently be observed. 


4.7.6. Non-local quantities of the body. A procedure statement written outside 
the scope of any non-local quantity of the procedure body is undefined. 


4.7.7. Parameter delimiters. All parameter delimiters are understood to be equi- 
valent. No correspondence between the parameter delimiters used in a procedure 
statement and those used in the procedure heading is expected beyond their 
number being the same. Thus the information conveyed by using the elaborate 
ones is entirely optional. 

4.7.8. Procedure body expressed in code. The restrictions imposed on a proce- 
dure statement calling a procedure having its body expressed in non-ALGOL code 
evidently can only be derived from the characteristics of the code used and the 
intent of the user and thus fall outside the scope of the reference language. 


5. Declarations 


Declarations serve to define certain properties of the identifiers of the pro- 
gram. A declaration for an identifier is valid for one block. Outside this block 
the particular identifier may be used for other purposes (cf. section 4.1.3). 

Dynamically this implies the following: at the time of an entry into a block 
(through the begin since the labels inside are local and therefore inaccessible from 
outside) all identifiers declared for the block assume the significance implied by 
the nature of the declarations given. If these identifiers had already been defined 
by other declarations outside they are for the time being given a new significance. 
Identifiers which are not declared for the block, on the other hand, retain their 
old meaning. 

At the time of an exit from a block (through end, or by a go to statement) 
all identifiers which are declared for the block lose their significance again. 

A declaration may be marked with the additional declarator own. This has 
the following effect: upon a reentry into the block, the values of own quantities 
will be unchanged from their values at the last exit, while the values of declared 
variables which are not marked as own are undefined. Apart from labels and 
formal parameters of procedure declarations and with the possible exception of 
those for standard functions (cf. sections 3.2.4 and 3.2.5) all identifiers of a 
program must be declared. No identifier may be declared more than once in 
any one block head. 
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Syntax. 
<declaration > :: = <type declaration)| {array declaration)| 
<switch declaration)| <procedure declaration) 


5.1. Type declarations 


5.1.4. Syntax. 
<type list>:: = <simple variable)| <simple variable), <type list) 
<type>:: = real | integer | Boolean 
<local or own type)>:: = <type>| own <type) 
<type declaration :: = <local or own type) <type list) 
5.1.2. Examples. integer #, g, s 
own Boolean Acryl, n 


5.1.3. Semantics. Type declarations serve to declare certain identifiers to re- 
present simple variables of a given type. Real declared variables may only 
assume positive or negative values including zero. Integer declared variables 
may only assume positive and negative integral values including zero. Boolean 
declared variables may only assume the values true and false. 

In arithmetic expressions any position which can be occupied by a real de- 
clared variable may be occupied by an integer declared variable. 

For the semantics of own, see the fourth paragraph of section 5 above. 


5.2. Array declarations 
5.2.4. Syntax. 


<lower bound): : = <arithmetic expression) 
<upper bound): : = <arithmetic expression 
<bound pair>:: = «lower bound): <upper bound) 
<bound pair list>:: = <bound pair}| <bound pair list), <bound pair) 
<array segment) :: = <array identifier) [<bound pair list >] | 
<array identifier), <array segment) 
<array list):: = <array segment)| <array list), (array segment) 
<array declaration): : = array <array list)| 
<local or own type> array <array list) 
5.2.2.. Examples. array a, b, c [7:n, 2:m], s[— 2:10] 
own integer array A[i/ c<0 then 2 else 1:20] 
real array g[— 7: —1] 


5.2.3. Semantics. An array declaration declares one or several identifiers to 
represent multidimensional arrays of subscripted variables and gives the dimen- 
sions of the arrays, the bounds of the subscripts and the types of the variables. 
5.2.3.1. Subscript bounds. The subscript bounds for any array are given in the 
first subscript bracket following the identifier of this array in the form of a bound 
pair list. Each item of this list gives the lower and upper bound of a subscript 
in the form of two arithmetic expressions separated by the delimiter: The bound 
pair list gives the bounds of all subscripts taken in order from left to right. 
5.2.3.2. Dimensions. The dimensions are given as the number of entries in the 
bound pair lists. ; 
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5.2.3.3. Types. All arrays declared in one declaration are of the same quoted 
type. If no type declarator is given the type real is understood. 


5.2.4. Lower upper bound expressions. 


5.2.4.1. The expressions will be evaluated in the same way as subscript expressions 
(cf. section 3.1.4.2). 

5.2.4.2. The expressions can only depend on variables and procedures which 
are non-local to the block for which the array declaration is valid. Consequently 
in the outermost block of a program only array declarations with constant 
bounds may be declared. 

5.2.4.3. An array is defined only when the values of all upper subscript bounds 
are not smaller than those of the corresponding lower bounds. 

5.2.4.4. The expressions will be evaluated once at each entrance into the block. 
5.2.5. The identity of subscripted variables. The identity of a subscripted 
variable is not related to the subscript bounds given in the array declaration. 
However, even if an array is declared own the values of the corresponding sub- 
~ scripted variables will, at any time, be defined only for those of these variables 
which have subscripts within the most recently calculated subscript bounds. 


5.3. Switch declarations 
5.3.14. Syntax. 
<switch list) :; = ¢designational expression) | 
<switch list), <designational expression) 
<switch declaration) :: = switch <switch identifier): = <switch list) 
5.3.2. Examples. switch S:=S1, S2, Q[m], if v>—5 then S32 else S4 

switch 0: =/1, w . 

5.3.3. Semantics. A switch declaration defines the values corresponding to a 
switch identifier. These values are given one by one as the values of the designa- 
tional expressions entered in the switch list. With each of these designational 
expressions there is associated a positive integer, 7, 2,..., obtained by counting 
the items in the list from left to right. The value of the switch designator corre- 
sponding to a given value of the subscript expression (cf. section 3.5. Designational 
expressions) is the value of the designational expression in the switch list having 
this given value as its associated integer. 
5.3.4. Evaluation of expressions in the switch list. An expression in the switch 
list will be evaluated every time the item of the list in which the expression 
occurs is referred to, using the current values of all variables involved. 
5.3.5. Influence of scopes. Any reference to the value of a switch designator from 
outside the scope of any quantity entering into the designational expression 
for this particular value is undefined. 


5.4. Procedure declarations 
5.4.4. Syntax. 
<formal parameter) :: = <identifier)> 
<formal parameter list)>:: = <formal parameter)| 
<formal parameter list) <parameter delimiter) <formal parameter) 
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<formal parameter part>:: = <empty)| (<formal parameter list) 

<identifier list) :: = identifier) | <identifier list), <identifier) 

<value part>:: = value <identifier list); | empty) 

<specifier>:: = string | <type>| array | <type) array | label | switch | 
procedure | <type> procedure 

<specification part >:: = <empty)| <specifier) {identifier list} ;| 
«specification part) (specifier) <identifier list; 

<procedure heading: : = <procedure identifier) <formal parameter part); 
<value part» <specification part) 

<procedure body): : = ¢statement)| (code) 

<procedure declaration>:: = procedure <procedure heading) (procedure body)| 
<type> procedure <procedure heading) «procedure body) 


5.4.2. Examples (see also the examples. at the end of the report). 
procedure Spur (a) Order: (n) Result: (s); value n; 

array a; integer m; real s; 

begin integer ; 

s:=0; 

for k: =1 step J until n do s:=s +a[k, k] 

end 


procedure Transpose (a) Order: (n); value n; 
array a; integer n; 
begin real w; integer 7, k; 
for 1: = 1 step J until n do 
for k: = 1+ step J until » do 
begin w: =a[it, k]; 
a(t, k]: =a[k, 2]; 
a[k,1|:=w 
end 
end Transpose 
integer procedure Step (u); real u; 
Step: = if OS uAuS/1 then / else 0 


procedure Absmax (a) size: (n,m) Result: (y) Subscripts: (i, k); 

comment The absolute greatest element of the matrix a, of size n by m is transferred 
to y, and the subscripts of this element to i and k; 

array a; integer n, m, 1, k; real y; 

begin integer #, @; 

y:=0; 

for p: =1 step J until » do for g: =1 step J until m do 

if abs(a[p, g])>y then begin y: =abs(a[p, q]); 1: =p; k: =q end end Absmax 
procedure Innerproduct (a, b) Order: (k, p) Result: (y); value k; 

integer k, p; real y, a, b; 


begin real s; 

$:=0; 

for p: =1 step J until k do s:=s+a x); 
y:=s 


end [nnerproduct 
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5.4.3. Semantics. A procedure declaration serves to define the procedure asso- 
ciated with a procedure identifier. The principal constituent of a procedure 
declaration is a statement or a piece of code, the procedure body, which through 
the use of procedure statements and/or function designators may be activated 
from other parts of the block in the head of which the procedure declaration 
appears. Associated with the body is a heading, which specifies certain identifiers 
occurring within the body to represent formal parameters. Formal parameters in 
th procedure body will, whenever the procedure is activated (cf. section 3.2.Func- 
tion designators and section 4.7. Procedure statements) be assigned the values 
of or replaced by actual parameters. Identifiers in the procedure body which 
are not formal will be either local or non-local to the body depending on whether 
they are declared within the body or not. Those of them which are non-local 
to the body may well be local to the block in the head of which the procedure 
declaration appears. 
5.4.4. Values of function designators. For a procedure declaration to define the 
value of a function designator there must, within the procedure body, occur an 
assignment of a value to the procedure identifier, and in addition the type of 
this value must be declared through the appearance of a type declarator as the 
very first symbol of the procedure declaration. 

Any other occurrence of the procedure identifier within the procedure body 
denotes activation of the procedure. 
5.4.5. Specifications. In the heading a specification part, giving information 
about the kinds and types of the formal parameters by means of an obvious 
notation, may be included. In this part no formal parameter may occur more 
than once and formal parameters called by name (cf. section 4.7.3.2) may be 
omitted altogether. 
5.4.6. Code as procedure body. It is understood that the procedure body may 
be expressed in non-ALGOL language. Since it is intended that the use of this 
feature should be entirely a question of hardware representation, no further rules 
concerning this code language can be given within the reference language. 


Examples of procedure declarations 
Example 1 
procedure euler (fct, sum, eps, tim); value eps, tim; integer tim; real procedure 
fct; real sum, eps; 
comment euler computes the sum of fct(t) for 1 from zero up to infinity by means of 
a suitably refined euler transformation. The summation is stopped as soon as tim 
times in succession the absolute value of the terms of the transformed series are found 
to be less than eps. Hence, one should provide a function fct with one integer argument, 
an upper bound eps, and an integer tim. The output is the sum sum. euler ts parti- 
cularly efficient in the case of a slowly convergent or divergent alternating series; 
begin integer i, k, n, t; array m[0:15]; real mn, mp, ds; 
t:=n:=t:=0;m(0]: =fct(0); sum: =m(0)/2 ; 
nextterm: 7: =1i+1; mn: = fct(t); 
for k: =0 step J until n do 
begin mp : =(mn-+ m[k])/2;m[k]:=mn;mn:=mp end means; 
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if (abs (mn) <abs(m[n])) A(n <15) then 
begin ds: =mn/2;n:=n+1;m[n]:=mn end accept 
else ds: =mn; 
sum :=sum+ds; 
if abs(ds)<eps then ¢:=it+1 else ¢: =0; 
if t<tim then go to nextterm 
end euler 


Example 21 
procedure RK(x,y,,F KT, eps, eta, xE, yE, fi); value x, y; integer 1; 
Boolean ji; real x, eps, eta, xE; array y, yE; procedure FKT; 
comment: RK integrates the system y,=f,(%, Vy, Vo, «++» Yn) (R=, 2, ... m) 
of differential equations with the method of Runge-Kutta with automatic search for 
appropriate length of integration step. Parameters are: The initial values x and y[k] 
for x and the unknown functions y,(x). The order n of the system. The procedure 
FKT(x, y, n, 2) which represents the system to be integrated, 1.e. the set of functions f,. 
The tolerance values eps and eta which govern the accuracy of the numerical integra- 
tion. The end of the integration interval xE. The output parameter yE which re- 
presents the solution at x =xE. The Boolean variable fi, which must always be 
given the value true for an isolated or first entry into RK. If however the functions y 
must be available at several meshpoints x, x,,..., %,, then the procedure must be 
called repeatedly (with x =x,, xE =x,,,, fork =O, 1,...,n—1) and then the later 
calls may occur with fi =false which saves computing time. The input parameters 
of FKT must be x, y, n, the output parameter z represents the set of derivatives z[k] = 
f,(*, v[1], y[2],..., y["]) for x and the actual y’s. A procedure comp enters as a 
non-local identifier ; 
begin 
array z, yl, y2, y3[1:n]; real x1, x2, x3, H; Boolean out ; 
integer k,7; own real s, Hs; 
procedure RKIST{(x, y,h, xe, ye); real x,h, xe; array y, ye; 
comment: RKIST integrates one single Runge-Kutta step with initial 
values x, y[k]| which yields the output parameters xe =x-+-h and ye[k], 
the latter being the solution at xe. 
Important: the parameters n, FKT, z enter RK1ST as non-local entities ; 
begin 
array w[1:n],a[1:5]; integer k,7; 
a[1):=a[2]:=a[5]:=A/2; a[3]:=a[4]:=h; xe: =x; 
for k: =1 step J until n do ye[k]: =w[k|:=y[k]; 
for 7: =I step J until 4 do 





1 This RK-program contains some new ideas which are related to ideas of S. GILL, 
A process for the step by step integration of differential equations in an automatic 
computing machine. Proc. Camb. Phil. Soc. Vol. 47 (1951) p. 96, and E. Fr6BERG, 
On the solution of ordinary differential equations with digital computing machines, 
Fysiograf. Sallsk. Lund, Férhd. 20 Nr. 11 (1950) p. 136—152. It must be clear how- 
ever that with respect to computing time and round-off errors it may not be optimal, 
nor has it actually been tested on a computer. 
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begin 
FKT (xe, w, n, 2); 
xe: =x+alj]; 
for k: =1 step J until » do 
begin 
w[k]:=y[k]+alj] x2[k]; 
ye[k] : =ye[h] +a[j+1] x2[k]/3 
end k 
end 7 
end RK1ST; 
Begin of program: 
if fi then begin H : = xE —x;s:=0Oend else H:=Hs; 
out : = false; 


AA: if (x+2.01xH —xE>0) =(H>0) then 
begin Hs: =H; owt : = true; H: =(xE —x)/2 end i; 
RKIST(x, y,2xH, x1, yl); 


BB: RKIST(x, y, H, x2, y2); RK1ST(x2, y2, H, x3, y3); 
for k: =1 step J until » do 
if comp (y1[k], y3[k], eta)>eps then go to CC; 
comment: comp(a, b,c) is a function designator, the value of which is the 
absolute value of the difference of the mantissae of a and b, after the exponents 
of these quantities have been made equal to the largest of the exponents of the 
originally given parameters a, b,c; 
x: =23; if out then go to DD; 
for k: =1 step J until n do y[k]:=y3[k]; 
if s=5 then begin s: =0; H: =2xH end 7}; 
s:=s+1; go to AA; 


CC: H:=0.5xH; out: = false; x1:=x2; 
for k: =1 step 1 until » do y1[k]:=y2[k]; 


go to BB; 
DD: for k: =1 step 1 until n do yE[k] : =y3[k] 
end RK 


Alphabetic index of definitions of concepts and syntactic units 
All references are given through section numbers. The references are given 
in three groups: 
def Fcllowing the abbreviation ‘‘def’’ reference to the syntactic definition 
(if any) is given. 


synt Following the abbreviation ‘‘synt” references to the occurrences in metalin- 
guistic formulae are given. References already quoted in the def-group are 
not repeated. 


text Following the word ‘‘text”’ the references to definitions given in the text 
are given. 
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The basic symbols represented by signs other than underlined (bold faced. 
Publishers remark) words have been collected at the beginning. The examples 
have been ignored in compiling the index. 
see: plus 
see: minus 
see: multiply 
— see: divide 
see: exponentiation 
Ss = 2 > = see: <relational operator) 

> VA — see: <logical operator) 
see: comma 
see: decimal point 
io see: ten 
: see: colon 
; see: semicolon 
:= see: colon equal 
ul see: space 
() see: parentheses 
[ ] see: subscript bracket 
‘’ see: string quote 
<actual parameter), def 3.2.1, 4.7.4 
<actual parameter list), def 3.2.1, 4.7.1 
<actual parameter part), def 3.2.1, 4.7.4 
<adding operator), def 3.3.4 
alphabet, text 2.1 
arithmetic, text 3.3.6 
<arithmetic expression, def 3.3.1 synt 3, 3.4.1, 3.3.1, 3.4.1, 4.2.1, 4.6.1, 5.2.4 
‘ text 3.3.3 
<arithmetic operator), def 2.3 text 3.3.4 
array, synt 2.3, 5.2.1, 5.4.4 
array, text 3.1.4.1 
<array declaration), def 5.2.1 synt 5 text 5.2.3 
<array identifier), def 3.1.1 synt 3.2.1, 4.7.1, 5.2.4 text 2.8 
<array. list), def 5.2.4 
<array segment), def 5.2.1 
<assignment statement), def 4.2.1 synt 4.1.1 text 1, 4.2.3 
<basic statement), def 4.1.1 synt 4.5.1 
<basic symbol), def 2 
begin, synt 2.3, 4.1.1 
<block), def 4.1.1 synt 4.5.4. text 1, 4.1.3, 5 
<block head), def 4.1.4 
Boolean, synt 2.3, 5.1.1 text 5.1.3 
<Boolean expression), def 3.4.1 synt 3, 3.3.1, 4.2.1, 4.5.1, 4.6.1 text 3.4.3 
<Boolean factor), def 3.4.1 
«Boolean primary), def 3.4.1 
<Boolean secondary), def 3.4.1 
<Boolean term), def 3.4.1 


WAzy~x 1+ 


_ 
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<bound pair), def 5.2.4 
<bound pair list), def 5.2.1 


<bracket), def 2.3 ‘ 

<code>, synt 5.4.1 text 4.7.8, 5.4.6 ‘ 
colon: , synt 2.3, 3.2.4, 4.1.4, 4.5.4, 4.6.4, 4.7.4, 5.2.4 | 
colon equal : = , synt 2.3, 4.2.1, 4.6.1, 5.3.4 


comma , , synt 2.3, 3.1.1, 3.2.1, 4.6.1, 4.7.1, 5.4.1, 5.2.4, 5.3.1, 5.4.4 
comment, synt 2.3 
comment convention, text 2.3 
<compound statement), def 4.1.1 synt 4.5.1 text 4 
<compound tail), def 4.1.1 
<conditional statement», def 4.5.1 synt 4.1.1 text 4.5.3 
<decimal fraction), def 2.5.4 
<decimal number), def 2.5.1 text 2.5.3 
decimal point ., synt 2.3, 2.5.1 
<declaration), def 5 synt 4.1.1 text 1, 5 (complete section) 
<declarator>, def 2.3 
<delimiter>, def 2.3 synt 2 
<designational expression), def 3.5.4 synt 3, 4.3.1, 5.3.4 text 3.5.3 
<digit>, def 2.2.1 synt 2, 2.4.1, 2.5.4 
dimension, text 5.2.3.2 
divide / —, synt 2.3, ¥%3.1 text 3.3.4.2 
do, synt 2.3, 4.6.1 
<dummy statement), def 4.4.1 synt 4.1.1 text 4.4.3 
else, synt 2.3, 3.3.1, 3.4.1, 3.5.1, 4.5.4 text 4.5.3.2 
<empty>, def 1.1 synt 2.6.1, 3.2.1, 4.4.1, 4.7.1, 5.4.1 
end, synt 2,3. 4.1.1 
entier, text 3.2.5 
exponentiation ¢ , synt 2.3, 3.3.1 text 3.3.4.3 
<exponent part), def 2.5.1 text 2.5.3 
<expression), def 3 synt 3.2.1, 4.7.1 text 3 (complete section) 
<factor)>, def 3.3.4 
false, synt 2.2.2 
for, synt 2.3, 4.6.1 
<for clause), def 4.6.1 text 4.6.3 
<for list), def 4.6.1 text 4.6.4 
<for list element), def 4.6.1 text 4.6.4.1, 4.6.4.2, 4.6.4.3 
<formal parameter), def 5.4.1 text 5.4.3 
<forma] parameter list), def 5.4.1 
<formal parameter part), def 5.4.4 
<for statement), def 4.6.1 synt 4.1.1, 4.5.1 text 4.6 (complete section) 
<function designator), def 3.2.1 synt 3.3.1, 3.4.1 text 3.2.3, 5.4.4 
go to, synt 2.3, 4.3.1 
<go to statement), def 4.3.1 synt 4.1.1 text 4.3.3 
<identifier), def 2.4.1 synt 3.1.1, 3.2.1, 3.5.1, 5.4.4 text 2.4.3 
<identifier list), def 5.4.1 
if, synt 2.3., 3.3.4, 4.5.4 
<if clause), def 3.3.1, 4.5.4 synt 3.4.1, 3.5.1 text 3.3.3, 4.5.3.2 
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<if statement), def 4.5.1 text 4.5.3.4 
<implication), def 3.4.1 
integer, synt 2.3., 5.1.1 text 5.1.3 
<integer», def 2.5.1 text 2.5.4 
label, synt 2.3, 5.4.1 
<label>, def 3.5.1 synt 4.1.1, 4.5.1, 4.6.1 text 1, 4.1.3 
<left part), def 4.2.4 
<left part list>, def 4.2.1 
<letter>, def 2.1 synt 2, 2.4.1, 3.2.1, 4.7.4 
<letter string), def 3.2.1, 4.7.1 
local, text 4.1.3 
<local or own type), def 5.1.1 synt 5.2.1 
<logical operator), def 2.3 synt 3.4.4 text 3.4.5 
<logical value), def 2.2.2 synt 2, 3.4.4 
<lower bound), def 5.2.1 text 5.2.4 
non-local, text 4.1.3 
minus —, synt 2.3, 2.5.1, 3.3.1 text 3.3.4.1 
multiply x, synt 2.3, 3.3.4 text 3.3.4.4 
<multiplying operator), def 3.3.4 
<number), def 2.5.1 text 2.5.3, 2.5.4 
<open string), def 2.6.1 
<operator>, def 2.3 
own, synt 2.3, 5.1.1 text 5, 5.2.5 
<parameter delimiter), def 3.2.1, 4.7.1 synt 5.4.1 text 4.7.7 
parentheses (), synt 2.3, 3.2.1, 3.3.1, 3.4.1, 3.5.1, 4.7.1, 5.4.4, text 3.3.5.2 
plus +, synt 2.3, 2.5.1, 3.3.4 text 3.3.4.1 
<primary>, def 3.3.1 
procedure, synt 2.3, 5.4.1 
<procedure body), def 5.4.1 
<procedure declaration), def 5.4.1 synt 5 text 5.4.3 
<procedure heading», def 5.4.1 text 5.4.3 
<procedure identifier), def 3.2.1 synt 3.2.1, 4.7.1, 5.4.1 text 4.7.5.4 
<procedure statement), def 4.7.1 synt 4.1.1. text 4.7.3 
program, text 1 
<proper string», def 2.6.1 
quantity, text 2.7 
real, synt 2.3, 5.1.1 text 5.1.3 
<relation>, def 3.4.1 text 3.4.5 
<relational operator), def 2.3, 3.4.1 
scope, text 2.7 
semicolon ; , synt 2.3, 4.1.1, 5.4.14 
<separator>, def 2.3 
<sequential operator), def 2.3 
<simple arithmetic expression), def 3.3.1 text 3.3.3 
<simple Boolean), def 3.4.4 
<simple designational expression, def 3.5.1 
<simple variable), def 3.1.1 synt 5.1.1 text 2.4.3 
space u, synt 2.3 text 2.3, 2.6.3 


Numer, Math. Bd. 2 10 
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«specification part), def 5.4.1 text 5.4.5 
<specificator)>, def 2.3 
<specifier>, def 5.4.1 
standard function, text 3.2.4, 3.2.5 
<statement>, def 4.1.1, synt 4.5.1, 4.6.1, 5.4.1 text 4 (complete section) 
statement bracket see: begin end 
step, synt 2.3, 4.6.1 text 4.6.4.2 
string, synt 2.3, 5.4.4 
<string>, def 2.6.1 synt 3.2.1, 4.7.1 text 2.6.3 
string quotes ‘’ , synt 2.3, 2.6.1, text 2.6.3 
subscript, text 3.1.4.1 
subscript bound, text 5.2.3.1 
subscript brackets [ ], synt 2.3, 3.1.1, 3.5.1, 5.2.4 
<subscripted variable), def 3.1.1 text 3.1.4.1 
<subscript expression), def 3.1.1 synt 3.5.1 
<subscript list>, def 3.1.4 
successor, text 4 
switch, synt 2.3, 5.3.4, 5.4.4 
<switch declaration), def 5.3.1 synt 5 text 5.3.3 
<switch designator), def 3.5.4 text 3.5.3 
<switch identifier), def 3.5.4 synt 3.2.1, 4.7.1, 5.3.4 
<switch list), def 5.3.1 
<term), def 3.3.1 
ten i, synt 2.3, 2.5.4 
then, synt 2.3, 3.3.1, 4.5.4 
transfer function, text 3.2.5 
true, synt 2.2.2 
<type>, def 5.1.4 synt 5.4.1 text 2.8 
<type declaration), def 5.1.1 synt 5 text 5.1.3 
<type list), def 5.1.1 
<unconditional statement), def 4.1.1, 4.5.4 
<unlabelled basic statement), def 4.1.1 
<unlabelled block», def 4.1.4 
<unlabelled compound), def 4.1.4 
<unsigned integer), def 2.5.1, 3.5.1 
<unsigned number), def 2.5.4 synt 3.3.1 
until, synt 2.3, 4.6.1 text 4.6.4.2 
<upper bound), def 5.2.1 text 5.2.4 
value, synt 2.3, 5.4.1 
value, text 2.8, 3.3.3 
<value part), def 5.4.1 text 4.7.3.1 
«variable», def 3.1.1 synt 3.3.1, 3.4.1, 4.2.1, 4.6.1, text 3.1.3 
«variable identifier), def 3.1.4 
while, synt 2.3, 4.6.1 text 4.6.4.3 


Regnecentralen, 
Copenhagen-Valby, Denmark 


(Received March 3, 1960) 
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Handbook for automatic computation 


Preparation of a handbook for automatic computation, in five or more volumes, is now 
under way for publication by Springer-Verlag. It will appear in the series, ,,Grundlehren 
der Mathematischen Wissenschaften“ (‘‘Yellow Series‘). Editors are 

F. L. BAUER, Mainz 

A. S. HouSEHOLDER, Oak Ridge 

F. W. J. OLVER, Teddington 

H. RuTISHAUSER, Zurich 

K. SAMELSON, Mainz 

R. SAvER, MUNICH 

E. STIEFEL, ZURICH 


The purpose of the handbook is to provide a collection of tested algorithms for mathematical 
computations of all sorts: the solution of finite and of functional equations, methods of 
approximating functions, the evaluation of special functions, etc. These algorithms are to 
be written in ALGOL, hence will be usable on any machine for which a suitable translator 
is available, and even without a translator can be used as a model for programming. It is 
evident that such a collection could have no general utility unless written in some common 
program language. The descriptive language will be English. 

As plans now stand, the organization of the series will be as follows: Volume 1a will contain 
a description of the use of ALGOL, and Volume 1b a description of the structure of trans- 
lators. These introductory volumes are the only ones that will not be made up primarily of 
actual algorithms. Volume 2 will be devoted to the solution of finite equations, linear and 
nonlinear, including the determination of characteristic values and vectors of matrices. 
Volume 3 will be on functional equations, especially. differential equations, ordinary and 
partial, and integral equations. Volume 4 is concerned with methods of approximation, and 
Volume § the evaluation of particular functions. It is possible that certain algorithms, such 
as those for solving inequalities, for mathematical programming, for statistical computations, 
and the like, that do not seem to fall naturally in any of these areas, may be reserved for 
a sixth volume. Each algorithm is to be accompanied by enough explanatory information 
to make it understandable, along with whatever information is available on speed, accuracy, 
range, or, more generally, for judging the effectiveness of the algorithm for a given type of 
problem. In any event, only pretested algorithms will be published. 

Before the appearance of the volumes themselves, the algorithms will be prepublished in a 
series of supplements to the journal, ,,Numerische Mathematik‘. This is partly to make 
generally available each algorithm at the earliest possible time. But in addition to this, it 
provides the possibility for including in the handbook itself additional information, and even 
corrections, that might come in from users. 

Contributions are earnestly solicited. For the present, at least, these must necessarily be in 
the form of actual algorithms, along with information as to the extent and mode of testing 
‘the algorithm, estimates of accuracy, and experience in using it. Untested algorithms will 
not necessarily be rejected ipso facto, but their use must necessarily await actual test. As 
algorithms are published, information relating to published algorithms also will be welcomed. 
Contributions may be sent to any of the editors named above. 
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